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Preface 


This book deals with the most essential elements of computer graphics, namely 
analytic geometry and programming. It explains how programmers can use plotters 
and other graphic devices without discussing in detail how these devices work, and 
in what types they are now available. I hope that the reader will appreciate such a 
device-independent approach. In any case we thus avoid the confusion that 
sometimes arises if hardware and software topics are intertwined. 

Much attention is paid to the graphical representation of three-dimensional 
objects. The first three chapters contain a few interesting programs such as one for 
B-spline curve fitting, but these chapters are also preparatory for the rest of the 
book. Chapter 4 gives the traditional transformations for wire-frame models. In 
contrast to this well-known subject, Chapter 5 introduces an efficient method for 
hidden-line elimination which, as far as I know, is new. Like all other algorithms in 
the book, the method is worked out in a complete program, called HIDLINPIX. 
Some applications of this program are given in Chapter 6. Since some of the 
programs may well be useful in practice, the book offers more than the ‘principles’ 
promised in its title. On the other hand, there are also some programs of a 
somewhat playful nature, which have no practical significance in themselves. I hope 
that in these programs the reader will find principles that also apply to more realistic 
problems. 

АП programs in the book are expressed in the C language. This might seem 
curious, since Pascal is more frequently used for such purposes. Having pro- 
grammed in a great many languages for more than a quarter of a century I would 
rate Pascal as very good but C as excellent. (Since I have written successful Dutch 
textbooks for them, I wish a long life to both languages!) It is difficult to prove that 
braces ( ) are at least as readable as the keywords begin or end, but the difference 
in length is obvious. I mention this because for some programs in this book only the 
compactness of the C language made it feasible to list them completely. 

Up to now I have delayed answering the difficult question for whom the book is 
intended. I do not know the curricula of universities and other institutions all over 
the world well enough to recommend the book for term X in course Y. However, it 
seems that for anyone who teaches computer graphics, at least some parts of the 
book will be useful. For example, simple solid objects can be drawn in perspective 
by preparing an input file for HIDLINPIX manually, so that even those readers who 
do not program could benefit from the book. There are exercises at the end of each 
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chapter but I am convinced that any teacher in this field can easily add problems to 
his or her personal taste. 

In the professional world many people are developing software packages for 
computer aided design. I hope that the book will also be useful in this area, and that 
it will be a modest contribution to the development of good products. 


L. Ammeraal 


СНАРТЕК 1 


Introduction 


1.1 A MOTIVATION FOR GRAPHIC PROGRAMMING 


Students will sometimes doubt the relevance of the things they have to learn and 
they will ask their teachers questions about them. If little time is available, a brief 
but unsatisfactory answer to such questions will simply consist of a reference to the 
examination requiring that knowledge. Fortunately, it is most unlikely for such 
questions to be asked if the subject matter is related to computer graphics. 
Compared with the usual listings of a line printer, graphic computer output is very 
attractive to look at, and even those who do not agree with this are convinced that 
computer graphics has useful applications. 

Far more satisfying than looking at graphic computer output is producing such 
results oneself. This do-it-yourself point of view applies to all arts, but in computer 
graphics we are in the unique position of having an extremely accurate and 
hard-working slave at our disposal. This slave can produce literally all kinds of 
drawn pictures, unless we are unable to instruct it how to do it. Unfortunately, the 
latter is not an exception but a rule. Most computer users are unable to make the 
computer draw the pictures they want, even if they spend a lot of money on 
software. We have to live with this unsatisfactory situation because it is impossible 
for the average computer-user to write all the software he needs. We have to buy 
software written by others. This does not imply that it should be unwise to train 
students in programming. We should remember that before software is available it 
must be written. It is highly improbable that all the software we need will soon be 
available, so in future, programming will be necessary. This applies particularly to 
graphic software, since many programs that produce tables will be replaced with 
those which have graphic output. 

Even those users who do not program but buy all the software they need will 
benefit from some knowledge of programming. Students should deal with concepts 
that will not be quickly outdated but will be valid forever. An algorithm to decide 
whether or not a surface in space hides a line segment is such a fundamental 
concept, whereas technical details of a commercially available program package are 
not. Of course, sooner or later the users of such a package have to study its details 
thoroughly, but by then they should already have some insight into graphic 
algorithms. Such insight is best acquired by studying the algorithms both theoreti- 
cally and practically. It is therefore recommended to program yourself, even if you 
think you will not be programming in the future. 


1.2 GRAPHIC PROGRAMMING AND THE С LANGUAGE 
Natural languages, though interesting in many respects, are inadequate to express 
algorithms. Not surprisingly, algorithms are best written in an algorithmic language. 
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An algorithmic language which is ‘understood’ by computers is called a high-level 
programming language, or, briefly, a programming language. We shall conform to 
this usage, but stipulate that the programming language be truly ‘high-level’: not 
only machines but also human beings should understand the language as efficiently 
and as easily as possible. We do not ask that the language be easily readable for 
everyone. Mathematical symbols are unreadable for most eight-year-old children 
but are almost indispensable for mathematicians and engineers. In much the same 
way the programming language we shall use will not be appreciated by everyone, 
but it has proved to be extremely useful for professional programmers. We shall use 
the C language, not only because it is widespread but also because of the qualities of 
the language itself. To program in C one must be accurate; in C, logical errors will 
lead to syntactic errors less often than in Pascal, for example. Such logical errors in 
a C program may lead to wrong results or messages with technical errors that are 
hard to understand. In other words, mistakes could have bad effects. We have to 
realize this if we write C programs. If we read them we want our programs to be 
easily readable, which is a completely different aspect. People who do not use this 
language themselves will sometimes think C programs too cryptic to be readable. A 
good deal of confusion exists about readability. If one line of a C program 
corresponds to ten lines of Basic it is unfair to complain that a single line of C text is 
not as readable as one of those ten lines in Basic. If we try to understand a short C 
program and presume that it must be simple because it is short, we shall be 
disappointed. This should be remembered on many occasions in this book; a very 
short C program might be non-trivial and even interesting! 

For a systematic study of the C language the reader is referred to Kernighan and 
Ritchie (1978). As to explaining remarks on the C language, we shall restrict 
ourselves to notations that are very specific for it. More facts about the C language 
can be found in the Appendix. Here is our first graphic C program: 


/* SQUARES: This program draws 50 squares inside each other */ 
main() 
{ float xà, yA, xB, yB, xC, yC, xD, yD, 
xxA, yyA, xxB, yyB, xxC, уус, xxD, уур, р, qi 
int 1; 
p=0. 95; 4=1. 0-р; 
xA=2. Oi хВ=8. 0; xC-8.0; xD=2. Oi 
уА=0. 5; ЧВ=0. 5; yC=6. 5; yD=6. 5; 
initgr(); 
for (1=0; 150; i++) 
€ move(xA, ЧА); 
drau(xB, yB); drau(xC, yC); drau(xD, yD); drau(xA, ЧА); 
xxA-p*xAá*q*xB; gygA-p*yA*q*uB; xxB=p#xB+q#xC; yyB=p#yB+q#uyC; 
xxC=p#xC+q#xD; yyC=p#yC+q#yDi xxD=p#xD+q#xAi yyD=p#yD+q#uyA; 
хА=ххА; xB=xxB; xC=xxCi xD=xxD; 
yA=yyAi yB=yyBi yC=yyCi yD=yyD 


endgr di 


The output of this program is shown in Fig. 1.1. 
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Fig. 1.1. Output of program SQUARES 


The program contains calls of four graphic routines: 


initgr( ) initializes graphic output; 

move(x, y) moves a (real or fictitious) pen to point (x, y); 

draw(x, y) draws a line segment from the current position of the pen to point 
(x, у); 

endgr( ) performs final actions such as clearing an output buffer. 


In C the technical term for ‘procedure’ or ‘(sub)routine’ is ‘function’. Function 
calls are written with parentheses, even if there are no arguments. The call initgr( ) 
is required before the functions move and draw can be called. Similarly, the final 
call of move or draw must be followed by a call of endgr. Both calls move(x, y) and 
draw(x, y) move a (real or fictitious) pen to point (x, y); with move(x, y) the pen is 
up and with draw(x, y) it is down. The four functions mentioned do not belong to 
the C language. They are external routines; after compilation of our program they 
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are added to it by a linking loader. If the routines are not directly available in this 
form they can easily be expressed in other graphic sub-routine calls that are 
available. 

We shall not be very specific about the hardware to be used. Instead of discussing 
various plotters and other graphic devices we shall focus on general, device 
independent programming principles. In the author's environment, the routine initgr 
asks the user to choose one of the options ‘Immediate’ or ‘Deferred’. The choice 
‘Immediate’ has the effect that graphic output is immediately sent to the screen of 
the graphic terminal, which offers the possibility of interaction. The option 
‘Deferred’ causes the plot information to be written onto a disk file, which 
afterwards can be converted and sent to several devices, such as various pen 
plotters, a matrix printer or the graphic terminal just mentioned. For all these 
devices everything works fine if the ranges for x and y are 0€ x <9, 0 « y <7. Thus 
point (x, y) is located x inches to the right of the y-axis and y inches above the 
x-axis, in other words, the origin O of the screen co-ordinate system is the 
bottom-left corner of the screen. The reader who is not happy with this will 
probably be able to apply a change of co-ordinates, especially after studying this 
book! Incidentally, in Chapter 2 the picture boundaries will be dealt with in a more 
elegant way. 

The graphical output shown in Fig. 1.1 consists of fifty squares. A square ABCD 
is drawn and then a new point A’ is chosen on side AB such that AA’ = 0.05 x AB. 
Points B', C' and D' are chosen similarly on the sides BC, CD and DA, 
respectively. The points A', B', C' and D' are now called A, B, C and D, and the 
procedure is repeated. 


EXERCISES 


1.1 Write a program to draw a great number of triangles inside each other similar 
to the way in which the squares are drawn in the program SQUARES. 
1.2 Write a program to draw the following line segments: 
From point (1.0, 6.0) to point (1.0, 1.0); 
From point (1.0, 5.8) to point (1.2, 1.0); 
From point (1.0, 5.6) to point (1.4, 1.0); 


From point (1.0, 1.0) to point (6.0, 1.0). 
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Two-dimensional algorithms 


2.4 TRANSFORMATIONS AND NEW CO-ORDINATES 
Consider the following system of equations 
х'=х+а 
* "=y 
We can interpret these equations in two ways: 


(1) All points in the xy-plane move a distance a to the right (see Fig. 2.1(a)). 
(2) The x- and y-axes move a distance a to the left (see Fig. 2.1(b)). 


This simple example shows a principle which also applies to more complex 
situations. We shall often deal with systems of equations, usually written as matrix 
products, and interpret them as a transformation of all points in a fixed co-ordinate 
system. However, the same system of equations can then be interpreted as a change 
of co-ordinates. 

We wish to rotate point P(x, y) through an angle ф about the origin O. The image 
point will be called P'(x', y') (see Fig. 2.2). Then there are numbers a, b, c, d such 
that x' and y' can be derived from x and y as follows: 

E = ax + by 


a.d 
у'=сх+ау (2.1) 


The values of a, b, с, d can be obtained by choosing first (x, у) = (1, 0). Setting 
x = 1 and y =0 in Eq. (2.1) we get: 
х'=а 
у'=с 
However, in this simple case the values of x' and y' are cos ф and sin g, as can be 
seen in Fig. 2.3(a). Thus we have: 
a = COS ф 
c=sing 
In the same way, Fig. 2.3(b) leads to: 
b = —sin ф 
d= coso 
Thus for Eqs (2.1) we now write: 


Poea (2.2) 


у' = x sin ф +усоѕ ф 
5 
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(a) (b) 


Fig. 2.1. (a) Translation; (b) change of co-ordinates 


0 x 


Fig. 2.2. Rotation about O through angle q 


(-sing, cos 9) (0, 1) 


(cos 9,sing) V 


0 (0, 1) x 0 х 
(а) (b) 
Fig. 2.3. (a) Image of (1, 0); (b) image of (0, 1) 


In the following program an arrow is rotated repeatedly through 6° about О and 
then drawn. 


/* QUADRANTI: #/ 
/# This program draws 14 arrows, flying counter-clockwise +/ 
/# about О in the first quadrant of the coordinate system  */ 


include <math. h> 
float x[41={ 6.0, 6.0, 9.9, 6.1), /* See Figure below #/ 
yl4]={-0.25, 0.25, 0.0, 0. 0O}; /* in this program */ 


TRANSFORMATIONS AND NEW CO-ORDINATES 


main() 


{ 


> 


int i; gi 


float pi, phi, cos phi, sin phi, xx, уу; 
рі=4. О+абап(1. 0); phi=6#pi/180; 
cos phi-cos(phi); sin_phi=sin(phi); 


initgr()i 
for (i=1i i<=14; i++) /* Initial position 
+ /* Rotate the arrow: #/ /* of the arrow: 
for (J=0i у<=3; j++) /* 
€ xx=x[jJ]; yy=y[j]; /* (6, 0.25) 
xLJl-xxs*cos phi-yyssin phi; ГАЈ 1 
yLjJ-xx*sin phi-*yys*cos phi; /* m 
+ /* FEN 
/* Draw the rotated arrow: +/ /+ 2/ _\3 
move(xtOl, yCO]); иж (5.9, 0): (6.1, О) 
for (y=1i Ј<=3; j++) /* i 
draw(xCyJ, y[j]); /+ О 
draw(x€1J, yC11); иж (6, -0. 25) 
} 
endgr(); 


Fig. 2.4. Output of program QUADRANT1 


*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
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In geometry points of an object are usually denoted by letters A, B,... Instead 
we use integers 0,1,... here. In the initial position the arrow points upwards; its 
centre is then located in (6,0). In the right-hand side of the program a 
not-too-successful attempt has been made to show the arrow in this position. The x- 
and y-co-ordinates of the ith vertex of an arrow are stored in the array elements x [i] 
and y[i] (i =0, 1, 2, 3). The arrays x and y are external: their definition is outside 
the function main. External arrays have the pleasing property that they can be 
initialized. Array subscripts count from 0 and we specify the number of elements. 
Hence the definition (also called declaration) 


float x[4] = (6.0, 6.0, 5.9, 6.1), ... 
introduces four array elements with their initial values: 
x[0]26.0  x[1]260 х[2]=5.9 х[3]=6.1 


Although initially y[0] is negative and all co-ordinates to be plotted must be 
positive, there is no problem, since, before the arrow is drawn, it is rotated, which 
brings it above the x-axis. Figure 2.4 shows the output of this program. 


2.2 ROTATION 


Equations (2.2) describe a rotation about O. Often this is not what we want. If we 
wish to rotate about a given point (xo, yo), we simply replace x with x — xo, y with 
y — yo, x’ with x' — x, and y’ with у’ — y, in these equations: 


x’ — xo = (x — xo) cos Фф — (y — Yo) sin y 
y' — Yo = (x — xo) sin ф + (y — yo) cos Фф 


pos cos ф— (y — yo) sin g (2.3) 


у' = yo + (x — xo) sin ф + (y — yo) cos p 


We apply this to our flying arrows: 


(+ ARROWS3O: */ 

/* This program draws 30 arrows, flying counter-clockwise #/ 

/# about point (xO, yO) */ 

#include <math. h> 

float x[41-4 0.0, 0.0, -0.08, O. 08}, /* See Figure belou */ 
yl4J={-0.25, 0.25, 0.0, 0.0 +; /* in this program #/ 

main() 

X int ir gi 


float pi, phi, cos phi, sin phi, dx, dy, 
хО=4. 5, yO=3. 5, r=3. Oi 
рі=4. Otatan(1.0);i рһі=12#рі/180; 
cos phi-cos(phi); sin phi-sin(phi); 
initgr(); 
/* Move to start position (xO+r, yO) : */ 
for (j20; JZ4; j++) € xLjl*-xO*r; yLj1*-gO; У 
for (і=0; 1<30; i++) 
€ /* Rotate the arrow: */ /* Initial position: */ 
for Су=0; J<£4i J++) (+ */ 
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€ dx2xLj1-xOi dy=ylyI-yOi ГАЈ (0, O. 12) */ 
xCy1=xO+dx#cos phi-dys*sin phi; /* 1 */ 
yCy1=yO+dx#sin_phi+dy#cos phi; /+# à #/ 
/# FEN */ 
/* Drau the rotated arrou: */ /* 2/ 3 #/ 
move(xCO], ytO1); /#(—0.04, O) i (0.04, O)#/ 
for (у=1; у<=З; j++) /# i */ 
drau(xL jJ, YC Jj); / i */ 
drau(xt11, yC11)i /# O #/ 
b /# (0, -0O. 12) */ 
endgr()i 
+ 
<a ګګ‎ 
AT FF. 


Fig. 2.5. Output of program ARROWS30 
Figure 2.5 shows the output of this program. 


2.3 MATRIX NOTATION 
Equations (2.2) can be written as one matrix equation: 


$ ow COS @ sing 
к yIB As pd (24) 
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x’ cos —sin X 
Мык Le es 

y sing соѕф Jly 
In books on computer graphics the row vector notation (2.4) is more frequently 
used than the column vector notation (2.5). We too shall use notation (2.4). In this 


notation the ith row of the square matrix is always the image of the ith unit vector 
(here i = 1, 2). It is possible to express Eqs (2.3) as a matrix equation: 


or, with column vectors: 


, 


cos sin 2 


2. 
—sing cos ф (2.6) 


k у]=[ »)*[x-xo y -»l| 
The right-hand side of this equation, however, is not a matrix product. In more 
complex situations, where rotations are combined with other transformations, it will 
be more convenient to have a single matrix product for each elementary transforma- 
tion. At first sight this seems impossible in cases where translations are involved. 
With the aid of 3 x 3 transformation matrices, however, it can be done, as we shall 
see. We begin with a simple translation. Point P(x, y) is translated to point 
P'(x', у’), where 


х'=х+а 
2,7 
ATTE (2.7) 
This could be written as 
1 0 
k у]=[х у 10 1 
a b 
but in view of what will follow we prefer the following equation: 
10 0 
к pt (2.8) 
x y' 1]=[x y Шо 1 0 
а b 1 


It is easily verified that Eqs (2.7) and (2.8) are equivalent. A technical term 
associated with this notation is ‘homogeneous co-ordinates’; this subject will be 
discussed in greater detail in Section 3.6. Writing every transformation as a 
multiplication of matrices will enable us to combine several transformations to one. 
To show such a concatenation of transformations we shall combine a rotation with 
two translations. A rotation through an angle ф about О was described by Eq. (2.4). 
We replace this equation with: 


cos9 sing 0 
[х у 1]=[х y 1]|-sing соѕф 0 
0 0 1 


(2.9) 


We shall now derive a new version of Eq. (2.6) to describe a rotation about 
(Xo, yo) through an angle ф; this equation will be of the form 


k »' l]J=[x y ЦА (2.10) 


where R is a 3 X 3 matrix. To find this matrix R we look upon the transformation as 
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a succession of the following three steps with (u,, vı) and (u2, v5) as intermediate 
points: 
(1) A translation to move (xo, yo) to О: 

(us vy l-[x y ЦГ 


where 
1 0 0 
Т'=| 0 1 0 
—Xo —Yo 1 


(2) A rotation through angle ф about О: 


[u> v; 1]=[u v; 1]Ro 
where 


cosp sing 0 Q.11) 
Ro= | -sing соѕф 0 ` 


0 0 1 
(3) A translation from O to (xo, yo): 
[x' y" i]J2[u v 1]Т 


where 
1 0 0 
Т=|0 1 0 
Xo yo 1 


The combination of these three steps is based on the fact that matrix multiplica- 
tion is associative, that is, 


(AB)C = A(BC) 


for any three matrices A, B and C whose dimensions are such that these 
multiplications are possible. For either side of this equation we simply write ABC. 
We now find: 


[x y’ 1]=[u v, ЦТ 
—-([u vi ЦАТ 
=[u, v, 1]RT 
={ у 1]T'}RoT 
= [х y 1]T'RT 
=[х y IR 
where 
R = T'RoT 
This is the desired matrix; performing two matrix multiplications gives: 
cosp sing 0 
К = | -sing соѕф 0 
Ci C> 1 
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where с; and c, are the following constants: 
C, = Xo — Xo COS P + yo sin ф 
C2 = Yo — Xo Sin 9 — yo cos ф 


2. WINDOWS AND VIEWPORTS 


In many situations we have to draw objects whose dimensions are given in units 
completely incompatible with our screen co-ordinate system. A building may be a 
hundred times larger than its image we want to produce. A molecule, on the other 
hand, is much smaller in reality than in a picture. Finally there are applications 
where the object to be drawn is not a concrete one but a graphical representation of 
relations between quantities, as, for example, in Fig. 2.6, where the profits of a 
certain company at the beginning of the twentieth century are shown. 

Problem-oriented dimensions are expressed in so-called world co-ordinates. In 
Fig. 2.6 the numbers 1901, 1902, 1903, 1904 and 50 000, 100 000, 150 000, 200 000 
and 250 000 are world co-ordinates. We now introduce the concept of a window; this 
is a rectangle surrounding the object (or a part of it) that we wish to draw, as shown 
in Fig. 2.6. The sides of a window are parallel to the x- and the y-axis. To avoid 
confusion in what follows, it is most important to note that a window is more closely 
related to the object than to the image to be produced. If, as usual, we introduce a 
horizontal x-axis and a vertical y-axis, the window in Fig. 2.6 is completely 
determined by: 


Xmin = 1898 
Xmax = 1908 

Ymin = — 150 000 
Ymax = 325 000 


We see that the dimensions and the position of the window are expressed in world 


Profit 


($) 
250 000 
200 000 
150 000 
] 100 000 
Window 
ee 50 000 


1901 1902 1903 1904 Year 


Fig. 2.6. Bar diagram in a window 
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co-ordinates. This may surprise the reader, since the window was introduced to 
specify what we wish to see in the picture and, at first sight, a number of inches 
seems more natural to achieve this than, for example, a fictitious profit of 
—$150000, as given here for умы. However, expressing windows in world co- 
ordinates is customary and very convenient in practice. 

To produce the desired picture a rectangular region of the screen must be given as 
well. This region is called a viewport. It is specified in a way similar to that of a 
window, that is, by the minimum and maximum values of the X- and Y-co- 
ordinates, but now they are screen co-ordinates. We shall denote them by capital 
letters X and Y. A typical viewport specification is: 


Xin = 1.5 
Xmax = 7.5 
Кш = 1.0 
Улы = 6.0 


The window will now be mapped to the viewport. For example, a given world 
co-ordinate x = 1898 will be converted to the screen co-ordinate X = 1.5. First, the 
scale factors f, and f, are calculated: 


f. Ex ax = X min 
" Xmax — Xmin 

f = T лах = Yin 
| Ymax — Ymin 


In our example we find f, = 0.6 and f, = 0.0000105. Then the distance X — X,,, of an 
image point to the left viewport edge is found as f, times the corresponding distance 
X — Хш Of the original point to the left window edge, and Y — Ymin is found 
similarly, which leads to: 


X- X min tf ` (х ارت‎ 


Y = Yain +, ” (y m i 


We conclude this section with three remarks. 


(1) The window may or may not encompass the complete object. If it does not, the 
parts of the object which are outside the window must not be drawn but they are 
to be cut off. This activity is known as clipping; it will be dealt with in Section 
2.5. 
In general the scale factors f, and f, are different. For a bar diagram this is just 
what we want, but it is not in cases where angles in the picture are to be the 
same as those in the object. We can then use the smaller of f, and f, for both 
scale factors. It is then recommended to replace Eq. (2.12) with formulae based 
on the centres of the window and the viewport. This will be shown in Section 
2.6. 
(3) The dimensions and the position of the window are not always known 
beforehand. In Section 2.6 they will be calculated instead of being specified by 
the user. 


(2 


МУ 
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2.5 LINE CLIPPING 


In this section we shall assume world co-ordinates and screen co-ordinates to be 
identical; in other words, the window and the viewport will coincide. The term 
‘window’ in this section could therefore be replaced with ‘viewport’ everywhere. It 
is, however, customary to clip against a window rather than against a viewport, and 
we shall conform to this usage. 

Figure 2.7 shows a rectangle ABCD which is a window. АП line segments to be 
drawn must be inside the window; in other words, if a line segment is to be drawn, 
those parts that are outside the window must be clipped. We wish the clipping 
process to be done automatically. Commands to draw triangle POR in Fig. 2.7 are 
to be interpreted as commands to draw the line segments P'P, PO and QQ'. 
Rectangle ABCD is drawn beforehand, so the complete polygon POQ'P' is drawn 
instead of triangle POR. 

, Since only the three points P, О and R are given, the co-ordinate pair (xp, уь) 
has to be derived from (xp, yp) and (xg, yp). In Fig. 2.7 we see that the slope of PR 
can be expressed in two ways, which gives the following equation: 


Ур’ TE YR — Yp 
Xp’ Е Хр Хв мары Хр 
We combine this with 


Ур' = Ymax 
and find 


(XR Хр)(Утах — yp) 

YR 7 Ур 
We see that the co-ordinates of P' can easily be calculated if it is known that 
endpoint P is inside the window and endpoint R(xp, yg) satisfies: 


Xp = Xp + 


Xmin < XR < Xmax 
JR > Ymax 


There are, however, many more cases to be considered. The logical decisions 


Fig. 2.7. Triangle to be clipped 
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0 Xmin Xmax x 


Fig. 2.8. Clipping in steps 


needed to find this out make line clipping an interesting topic from an algorithmic 
point of view. In Fig. 2.8 it is clearly not sufficient to clip line segment PQ against 
line CD. Cohen and Sutherland developed an algorithm for line clipping, which is 
presented in Pascal by Newman and Sproull (1979). We shall express this algorithm 
in the C language. 

With any point P(x, y) we associate a four-bit code 


bı b, b, bo 


where b, is either 0 or 1 (i = 0, 1, 2, 3). This code contains useful information about 
the position of P relative to window ABCD. In C the (truth) values of expressions as 
x < Xmin are 1 for true and 0 for false. Using this, we write 


b3 = (x < xmin) /* P to the left of AD  */; 
b2 = (x > xmax) I* P to the right of BC */; 
b1 = (y € ymin) /x P below AB ж/; 
Ь0 = (у > утах) /* Р above CD */; 
Only nine out of the sixteen possible bit-configurations actually occur, and these are 


shown in Fig. 2.9. 
In C such a code value is delivered by the following function: 


int code(x, y) float x, y; 
(return(x < xmin) «3 | (x > xmax) «2 | (y < ymin) «1 | (v > утах); 


) 


To understand this, we must know that b <n is what we obtain if the bitstring b is 
shifted n positions to the left. Besides << for left shift there is a bit operator written 
as |, for bitwise or. This operator must not be confused with the logical or-operator 
written as ||. A less efficient version of the above expression in the return statement 
would be: 


(x < xmin)*8 + (x > xmax)*4 + (y < ymin)*2 + (y > ymax) 
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Fig. 2.9. Code values 


The latter version is merely given to understand the former. 

The above function code will be used in the function clip, whose task is to analyse 
a given line segment P,P, and to draw that part of this line segment that is in 
window ABCD, if such a part exists. It works as follows. 

As long as at least one of the codes for P, and P, contains a 1-bit, either P, or P; is 
moved from outside the window to one of its edges or to an extension of such an 
edge. In the latter case the point is still outside the window, so another move will be 
necessary. In Fig. 2.8, for example, a move from P to R has to be followed by one 
from R to S. Then at the other end of the line segment clipping may be necessary, as 
Fig. 2.8 shows. Thus clipping is a repetitive process; in each step the distance 
between P, and P; decreases. The process terminates as soon as both points are no 
longer outside the window. Line segment P,P, thus obtained is then drawn. There 
is, however, another important case in which the loop is to terminate, namely if 
both P, and P; are outside the window and at the same side of the window. This 
case may not apply initially but may arise during the clipping process. If the 
endpoints of a line segment are outside a window, the line segments may or may not 
intersect the window, as Figs 2.8 and 2.10 show. In Fig. 2.10, initially P; and Р, are 
not both below the window. After some steps in the clipping process О and S are 


Fig. 2.10. Line outside window 
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the new positions of Р, and P}, respectively. Since both points are now below the 
window it can be decided that nothing has to be drawn at all. Such decisions are 
made on the basis of the value of code(x1, y1) and code(x2, y2). The points P, and 
P, are at the same side of the window if and only if their codes have a 1 in the same 
position. For the three points Р,, Q, S the third bit from the left (b1) in their codes 
is 1, whereas this bit is 0 for point P}. This is why P, has to be moved to S. Moving 
P, to О is not necessary but this move is performed because it does no harm and 
keeps the algorithm simple. 

Like the bitwise or-operator |, explained earlier, the C language offers the bitwise 
and-operator, written as &. Note that the logical versions of these operators are || 
and &&. Bit operators (such as & and |) give bitstrings as results; a logical operation 
gives a single 1 or 0, being the representation of true and false, respectively. Not 
only 1 but any value unequal to zero acts as true if it is used in a logical context. 
Therefore the loop construction 


while (c1|c2)... 


will execute the actions indicated by..., as long as evaluation of c1 |c2 gives a 
bitstring containing a 1-bit, that is, as long as there is a 1-bit anywhere in c1 or c2. 
On the other hand, the statement 


if (c1 & c2) return; 


causes a direct return from the function, if and only if the bitstring value of c1 & c2 
contains a 1-bit, that is, if and only if both cl and c2 have a 1-bit in the same 
position. The complete function clip is shown at the bottom of the following 
program, which clips concentric pentagons. 


/* CLIPDEMO: #/ 
/# Demonstration of the Cohen & Sutherland line-clipping #/ 
/# algorithm */ 


include <math. hè 
float xmin=1.0, xmax=7.0, ymin=2. 0, ymax=6. О; 


main() 
f int ii 
float r, pi, alpha, phiO, phi, xO, yO. х1, yl. х2, ya; 
pi=4. O*atan(1.0); alpha=72. О#рі/180. Oi phiO-O.O; 
x0-4.0; yO-4.0; 
initgr(); 
/* The window is nou draun: */ 
move(xmin, ymin); draw(xmax, ymin); drau(xmax, ymax); 
draw(xmin, ymax); draw(xmin, ymin); 
/* As far as permitted by the boundaries of the window, +#/ 
/# 20 concentric regular pentagons are drawn: */ 
for (r=0. 5; r<10. 5; т+=0.5) 
< x2=x0+r#cos(phiO); ye2=yO+r#sin(phiO); 
for (1=1; i<=5; i++) 
€ phi=phiO+i#alpha; 
х1=х2; y1=y2; 
x2=xO+r#cos(phi)i y2-yO*rssin(phi); 
clip(x1, уі, x2; y2); 
} 
} 
endgr()i 
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int code(x, у) float х, yi 
€ return (х<хтіп) <<3 i (х>хтах)<<2 i (у<утіп) <<1 | Cy>ymax)i 
T 


clip(x1, уі, хә, уә) float xl, уі, хе, ya; 
< int ci=code(xi, yl), c2-code(x2, уо); float dx, dy; 
while (ciic2) 
€ if (ci&c2) returni 
dx=x2-x1i dy=y2-yli 
if Cer? 


€ if (х1<хтїп) < yl += dys(xmin-x1)0/dxi xi-xmin; } else 
if (xi>xmax) X yl += dy*(xmax-x1)/dx; xl=xmaxi >} else 
if (yl<ymin) € x1 += dx#(ymin-y1)/dy; yl=ymini } else 
if (yi»ymax) < x1 += dx*(ymax-yl)/dyi; yl=ymaxi } 
ci=code(x1, у1); 

) else 

< if (x2€xmin) < ye += dy*#(xmin-x2)/dxi x2=xmini } else 
if (x2>xmax) < y2 += dy*(xmax-x2)/dx; x2=xmaxi } else 
if (y2@<ymin) < x2 += dx¥(ymin-y2)/dyi y2=ymini } else 
if (y2>ym- ) < x2 += dx*(ymax-y2)/dyi y2=ymaxi } 


ce=code(x2, y2); 
} 
У 
томе (хі, yl); draw(x2, ya); 


Figure 2.11 shows the output of this program. 


2 


Fig. 2.11. Output of program CLIPDEMO 
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2.6 ADJUSTING SIZE AND POSITION AUTOMATICALLY 


To have a picture drawn within the boundaries of a viewport we can first clip it 
against a given window, as in Section 2.5, and then map the window and its contents 
to the viewport, as in Section 2.4. For many applications this procedure is 
satisfactory. In this section our approach will differ from this in the following 
respects: 


(1) The object will be completely drawn, so no clipping will take place. 

(2) The window will be calculated rather than specified. 

(3) In mapping from window to viewport the same scale factor will be applied to the 
horizontal and the vertical directions. 


Point (1) requires that the object be finite. For most applications this will be no 
severe restriction but it will exclude landscapes. Point (2) can be realized by 
scanning the object data twice. During the first scan the window boundaries X mins 
Xmax» Ymin aNd Ymax are determined. The drawing will be produced during the second 
scan. We shall use a file on disk for this purpose to avoid trouble due to memory 
limitations. Point (3) implies that each triangle in the picture will be similar to its 
original triangle in the object, which means that the mapping leaves angles 
unaltered. 

Suppose that we are given the triangle POR of Fig. 2.12. The relevant world 
co-ordinates are: 


Xp=1.0  xg-1.5 Xg — 1.2 

Ур = 0.8 yo = 0.9 Ув =1.1 

This will lead to a calculated window with characteristics: 
Xmin = 1.0 Xmax= 1.5 
Утіп = 0.8 Ymax = 1.1 


Notice that extreme object points are now on the window boundaries, which was not 
the case in Section 2.4. To have some blank space at the four sides of our screen or 


R(1.2, 1.1) 


1 
A 0.9) 


Р(1.0, 0.8) 


1 2 


Fig. 2.12. An object to be adjusted 
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paper we shall have to specify a viewport that is somewhat smaller than it would 
otherwise have been. For example, we can set Xpin=0.2 instead of 0.0. Our 
viewport is completely determined if we choose, for example: 
Xmin=0.2  Xmax= 8.2 
Yin = 9.5 Ps = 6:9 
As in Section 2.4, we calculate: 
Lx Au. S2 =0 _ 
Xma — Xmin 15—10 
Fass las. 0-2—0Ә 
Í = = = 
Ymax ^ Ymin 1.1 — 0.8 
The smaller of f, and f, will now be our scale factor f. Recall that distances will be 
multiplied by the scale factor, so part of the picture would certainly have been 


outside the viewport if we had chosen a scale factor greater than f, or fy. In our 
example we have: 


Lh 16 


20 


f=f,=16 
It will be clear that this common scale factor causes the mapped triangle to fit 
exactly in the viewport in the x-direction but not so for the y-direction, where blank 
space is to be added. We wish to distribute this blank space equally between the 
lower and the upper part of the viewport. This is accomplished by using the central 
value Үс instead of the minimum value Ymin, as in Section 2.4, to calculate the 
constant c>. The value of c, is determined similarly: 


хс = 0.5(Xmin + Xmax) = 0.5(1.0 + 1.5) = 1.25 
ус = 0.5(ут + Ymax) = 0.5(0.8 + 1.1) = 0.95 


Xe =0.5(Xmin + Xmax) = 0.5(0.2 + 8.2) = 4.2 
Ye = 0.5(Y.i + Ymax) = 0.5(0.5 + 6.5) = 3.5 
e = Xc- f: xc242—16x125- -15.8 
с»= Үс f : yc = 3.5 - 16x 0.95 = 11.7 


For any point (x, y) of the object its image point (X, Y) is now calculated as: 


X=f-x+c,=16x — 15.8 
Y-f:y*c5716y – 11.7 


We shall now develop a program which draws a picture whose window cannot be 
specified beforehand and a number of new aspects of the C language will appear in 
this example. We shall use random numbers to generate a curve of unpredictable 
shape and size and, as usual, the curve will be approximated by a great many line 
segments. Automatic scaling and positioning will relieve us from a practically 
impossible task. We shall generate co-ordinates x and y for each line segment and 
write them to a file on disk. More precisely, we write so-called structures, containing 
the triples 


x y code 
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where (x, y) is the point a pen has to go to, and code is either 0 or 1, meaning ‘pen 
up’ or ‘pen down’, respectively. In other words, 


x у 0 means move(X, Y) 
x у 1 means draw(X, Y) 


where X and Y are screen co-ordinates corresponding to the world co-ordinates x 
and y. In this example we shall use a special program which generates a curve and 
writes triples in the file A.SCRATCH. The execution of this program is to be 
followed by a run of the general plot program GENPLOT, which reads the triples 
twice; first, to determine the window characteristics xmin, xmax, ymin, ymax and 
second, to perform the actual move and draw operations, using screen co-ordinates 
X and Y, derived from the world co-ordinates x and y. 

In the curve-generation program we start in the origin O and move one unit of 
distance at a time. There is always a current direction q and a current turning angle 
а. Initially they are both set to zero. Before each step, а is increased by a randomly 
chosen angle between —6° and +6° (integer). The new turning angle « is then added 
to the current direction ф to find a new direction. We shall limit the curvature and 
reduce the chance of looping in circles. Therefore the above algorithm is modified in 
that we set œ equal to zero as soon as its absolute value is greater than 15°. The 
following program generates the curve. 


/* CURVGEN: Generation of a random curve */ 
include <math. h> 
main() 
< int i, №500; 
float x-0.0, y=0.0, xO, yO. phi, direction(); 
pfopen(); pmove(x, у); 
for (i=l; i<=Ni i++) 
< xO=xi yO=y; 
phi-direction(); 
x-xO*cos(phi); y=yO+sin(phidi 
pdrau(x, у); 
} 
pfclose(); 
} 


float direction() 
X static int рһі=0, alpha=O, first=i; 
/# Static variables are initialized in the first call only! #/ 
float pi-3.1415926; long int seed; 
if (first) < first=0; time(£&seed); srand((int)seed); } 
alpha*-rand()X13-6; 
if (abs(alpha)>15) alpha=0; 
phi+=alpha; 
return ((float)phi#pi/180. О); 
+ 


#include <stdio.h> 

FILE #fpi 

struct {float xx; float yyi int iii} s; 
pfopen() < fp=fopen("a. scratch", "u"); } 


pmove(x, у) float x, yi 
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X s.xx-xi в.уу=у; s.ii=0; /* О = pen up #/ 
furite(&s, sizeof s, 1, fp); 
} 


pdrau(x, у) float x, yi 
{ s.xx-xi в.уу=у; s.ii-l; /* 1 
furite(&s, sizeof s, 1, #р); 


pen down #/ 
} 


pfclose() < fclose(fp); } 


The function direction shows a call of the random initialization function srand. Its 
argument seed supplies a start value for random number generation. The function 
time is used to obtain a value for seed depending on the actual clock time. In this 
way we generate different curves if the program is run more than once. The function 
rand gives a large non-negative integer. This is converted to one of the integers 
0,1,...,12 by taking the remainder after dividing it by 13. Thus 


O<rand( )%13 <12 
—6<rand( )7013— 6x6 


More new C aspects concern input and output (I/O). Most I/O functions require 
the line #include (stdio.h), which ‘includes’ a so-called header file for standard 
I/O. In C we distinguish formatted and unformatted I/O. The following functions 
are used very often: 


scanf: formatted input from the terminal, 
printf: formatted output to the terminal, 
fscanf: formatted input from a file, 
fprintf: formatted output to a file, 

fread: unformatted input (from a file), 
fwrite: unformatted output (to a file). 


Formatted I/O deals with readable characters; there is a line structure in the same 
way as on a printed page. Unformatted data on a file has the same structure as it has 
in memory; for example, an integer is represented by a fixed number of bits. We 
have used unformatted I/O for reasons of efficiency. A file is ‘opened’ by a call of 


29099 


the function fopen. Its second argument is either the string "7" to initiate reading or 
the string "w" for writing. (Readers who use a Prime computer must replace "r" 
with "i" and "w" with "o" for unformatted I/O.) Variable s is the structure 
containing the three numbers x, y and code that are to be written. A ‘pointer’ &s to 
this variable is the first argument of fwrite. We can regard &s as a notation for the 
address of s. The second argument sizeof s is the size of one structure s and the third 
argument 1 is the number or structures to be written. The fourth argument fp is the 
file pointer, obtained by the call of fopen. 

We now turn to the general program which reads triples, determines the window, 
performs the conversion from world co-ordinates to screen co-ordinates and finally 
draws the picture within a given viewport. The corners of the viewport are drawn 
too, along with a dot in the middle of the bottom boundary so that in abstract 
pictures we can distinguish between bottom and top. (If the dot and the corners are 
not desired they can easily be erased.) 
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Notice that in this program the file A.SCRATCH is opened twice. Closing and 
opening has the effect of rewinding a file: 


/# GENPLUT: A general adjusting and plotting program. */ 
/# The file A. SCRATCH contains input data. #/ 
#include <stdio.h> 


main() 
€ float x, y, xmin, xmax, ymin, ymax, X, Y, Xmin, Xmax, 
Ymin, Ymax: fx, fy, f, xC, ЧС, ХС, YC, cl, с2; 


FILE #fpi 
struct {float xx; float уу; int iii} si 
fp=fopen("a. scratch", "r"); /* "i" on PRIME computers 3/ 


xmin=ymin=1e30i xmax=ymax=-xmin; 
while (fread(&s, sizeof s, 1, fp)) 
€ x=5. xx; у=5. Цу; 
if (xCxmin) xmin=x; 
if (x>xmax) xmax-x; 
if Cy<ymin) ymin=yi 
if (y>ymax) ymax=yi 
> 
fclose(fp); 
init vieuport(£Xmin, &Xmax, &Ymin, %Үтах); 
fx=(Xmax-Xmin)/(xmax-xmin); fy=(Ymax-Ymin)/(ymax-ymin); 
f2(fxCfy?fx: FU); 
xC=0. S* xmintxmax); yC=0. 5#(ymin+ymax ); 
ХС=0. 5#(Xmin+Xmax); YC=0. 5#(Ymin+Ymax ); 
с1=ХС-##хС; c2=YC-f#yC; 
fp=fopen("a. scratch", "r"); /# "i" on PRIME computers #/ 
while (fread(&s, sizeof s, 1, fp)) 
€ x-s. хх; у=5. уу; 
Х=##х+с1; Y=f#y+c2; 
if (5.11) drau(X, YO; else move(X, Y); 
+ 
fclose(fp); endgr(); 
F 


init vieuport(pXmin, pXmax, pYmin, pYmax) 
float *pXmin, *pXmax, *pYmin, #pYmax; 

+ float Xmin, Xmax, Ymin, Ymax, ерѕ=0. 2; 
printf("Give viewport boundaries Xmin, Xmax, Ymin, Ymax\n"); 
scanf ("Af Xf Xf Xf", £Xmin, £Xmax, £Ymin, £Ymax)i; 


/* Shou the four vieuport corners: */ 

initgr(); 

move(Xmin, Ymin*eps); drau(Xmin, Ymin); drau(Xmin*eps, Ymin); 
move(Xmax-eps, Ymin); drau(Xmax, Ymin); drau(Xmax, Ymin*eps); 
move(Xmax, Ymax-eps); drau(Xmax, Ymax); drau(Xmax-eps, Ymax); 


move(Xmin*eps, Ymax); drau(Xmin, Ymax); draw(Xmin, Ymax-eps?); 
move((Xmin+Xmax)/2, Ymin); drau((Xmin*Xmax)/2, Ymin); 

/* Dot in the middle of bottom boundary for orientation #/ 
*#pXmin=Xmini 3pXmax-Xmax; #pY¥min=Ymini #pYmax=Yma x; 


The arguments of the function init.viewport are pointers. The notation float 
*pXmin expresses that *pXmin is of type float, so pXmin is of type pointer to float. 
The unary operators & and * are each other's inverses, so pXmin = &Хтїп and 
Xmin = *pXmin. The output of this program is shown in Fig. 2.13. 
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Fig. 2.13. Result of programs CURVGEN and GENPLOT 


2.7 APPLICATIONS OF RECURSION 


Many tasks can be formulated recursively. An example of such a task is: 
With three given numbers xc, yc and r connect the points with co-ordinates 


Xj = Xc +r cos Q; 

yi = yc + r sin Q; 

(i= 0, 1, 2, 3, 4, 5; p; =i · 144°) 
which will produce a star. Subsequently (and as part of the task!) perform a 
similar task five times, but now with the three given numbers: 

Xe: = хс + 2r cos ф; 

yc: = yc + 2r sin Q; 
r'=0.5r 
(/= 0, 1, 2, 3, 4; ф,= 36° +]. 72°) 


The whole task is to be performed only if the given value r is at least 0.1. 
For the initial or main task the three given numbers are xc=0, ус= 0, r=1. As 
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before, we shall use automatic size adjusting and positioning. Here is a program 
which performs the task described above: 


/* STARS: Stars of various sizes */ 
#include <math. h> 


main() 
{ pfopen(); star(O.O0, О. О, 1.0); pfclose(); 
} 


star(xC, ЧС, т) float xC: yC, Fi 
€ float phi, r half, rm double, 
factor=0. 0174533; /* factor = pi/180 #/ 
int ii 
if (reo. 1) returni 
pmove(xC+r, ЧС); 
for (1=1; iz-25; i++) 
€ phi=i#144#factor; 
pdraw(xC+r#cos(phid, yC+r#sin(phi)); 
F 
r_half=0. S#ri r_double=2#ri 
for (i20; 15; i++) 
< phi=(36+i#72)#factor; 
star(xC«*r doublescos(phi), уС+т doublessin(phi), r half); 
} 


#include <stdio.h> 
FILE #fp; 
struct {float xxi float yy; int iii} si 


pfopen() < fp=fopen("a. scratch", "ш"); } 


pmovetx, у) float х, yi 

€ s.xx-xi 5. yy=yi s.ii-O; /# О = pen up #/ 
furite(%s, sizeof s, 1, fp); 

} 


pdrau(x, y) float x, yi 

€ s.xx-xi S.yy=yi s.ii-l; /* 1 = pen down #/ 
furite(£s, sizeof s, 1, #р); 

Y 


pfclose() < fclose(fp); > 


This program, followed by the program GENPLOT of the previous section, will 
produce Fig. 2.14. 

Our next example is known as Pythagoras's Tree. This tree is often shown as in 
Fig. 2.15. Each right-angled triangle in this tree has an angle of 45°. We shall again 
use random numbers in a rather general program: it can produce the tree of Fig. 
2.15, but also less regular trees. The angles that are 45? in Fig. 2.15 will be randomly 
chosen between (45 — delta)? and (45 + delta)", where delta is given as input data, 
together with n, the recursion depth. The regular version of Fig. 2.15 was produced 
by choosing delta — 0 and n — 7. In the picture n is the number of triangles we 
encounter if we follow a path from the root to a leaf of the tree. The heart of the 
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Fig. 2.15. Pythagoras's Tree, regular version 
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Fig. 2.16. Point numbering 


program will be the recursive function square.and triangle, with the recursion depth 
n as its first argument. If n is greater than zero, the task of square and triangle is, as 
its name indicates, to draw a square and a triangle on top of it, and subsequently 
invoke the same function twice with appropriate new arguments, the first of which is 
n —1. The size and position of the square are fully determined by the four 
arguments Xo, yo, a and g (see Fig. 2.16). To draw the triangle we need the angle a. 
This angle, expressed in degrees, is equal to 45 + deviation, where deviation is one 
of the integers —delta, —delta + 1, . . . , delta, randomly chosen as in Section 2.6. 

In Fig. 2.16 the relevant points are numbered 0, 1, 2, 3, 4. The co-ordinates x, 
and y, of point 0 are given. To find the co-ordinates of the other points we first 
consider the much simpler situation with Фф = 0, that is, where side 0 1 of the square 
is horizontal. In this situation, the co-ordinates of the points are easily found. They 
are stored in the arrays x and y. Then we rotate everything about point 0 through 
the angle @ in the same way as in Section 2.3. The results of the rotation are stored 
in the arrays xx and yy: 


/# PYTH TREE: Variants of the tree of Pythagoras */ 
include <math. h> 

*define рі З. 1415927 

int deltai 

long int seed; 


main() 

< int ni 
pfopen(); time(&seed); srand((int)seed); 
printf("Give angle delta in degrees (О < delta < 45)\n"); 
scanf("Xd", &delta); 


printf("Give recursion depth n\n"); scanf("Zd", &n); 
square and triangle(n, О. О, О.О, 1.0, 0.0); 
pfclose(); 
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Square and triangle(n, xO, yO, a, phi) 
int n; float xO, yO, а, phii 
< float xt51, yL51, xxt51, yyL51, cphi, sphi 
alpha, calpha, salphai 
int 1, deviationi /* phi and alpha 
if (п==0) returni 


deviation=rand()%(2#delta+i)-delta; 
alpha=(45+deviation)#pi/180. 0; 
xCO]=xC3]=x0i xt11-2xt21-2xO-*ai 
yCOJ=yliI=yO; yt21-yt31-yO*a; 
calpha=cos(alpha)i salpha=sin(alpha)i 
c=atcalphai b=axsalpha; 
xC4J=xC3J+c#calphai 
yC4]1=yl3]1+c#salpha; 
/# Rotation about (xO, yO) 
/* this uas explained in Section 2-3 
cphi=cos(phi)i sphi=sin(phidi 
ci=xO-xO#cphi+yO#sphi; 
c2=y0-x0O#sphi-yO#cphi; 
for (1=0; 1<5; i++) 
€ xxCi]=xCil#cphi-ylil#sphi+ci; 
yyCil=xCil#sphi+ylil#cphi+ca; 


> 

ртоме(ххС3], yy[31); 

for (i20; 1<5; i++) pdrau(xxLilJ, yyliddi 
pdraw(xxf2], yyt21); 
square_and_triangle(n-1, xxC€3], yyl31, с, 
square_and_triangle(n-1, xxt42J, уу[41, b, 


} 


#include £stdio.h> 


FILE *fpi 


struct {float xxi float yyi int iii} si 

pfopen() < fp=fopen("a. scratch", "ш"); + 

pmove(x, у) float х, yi 

{ s. xxzxi 5.уу=у; 5. іі=0; /+# О = pen up */ 
furite(%s, sizeof s, 1, #р); 

} 

pdrau(x, у) float х, yi 

{ s. xxzxi 5. yy=yi 5. 11=1; /* 1 = pen down */ 
furite(%s, sizeof s, 1, #р); 

} 

pfclose() {£ fclose(fp); ) 


/* delta in degrees 


through angle phi; 
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cil, сә, b, с, 


in radians */ 
*/ 


*/ 
*/ 


phitalpha)i 
phi*alpha-O. S#pidi 


Again, the file A.SCRATCH, produced by this program, has to be processed by 
program GENPLOT of Section 2.5. The graphical output thus produced, with 


delta = 30 and n — 7, is shown in Fig. 2.17. 


2.8 CURVE FITTING 


In computer aided design (CAD) and computer aided manufacturing (CAM) it is 
often required to construct a smooth curve or a smooth surface through some given 
points. Here we shall deal with two dimensions only, so we shall restrict ourselves to 
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Fig. 2.17. Typical output of program PYTH.TREE 


curves in the xy-plane. Curve fitting will be a sound base for surface fitting at a later 
stage. 

Out of several methods available we choose the B-spline form. A sequence of 
points is given, and between two successive points of the sequence a cubic curve is 
constructed, based on the position of four points, namely the two just mentioned 
and their two neighbouring points. The B-spline form yields smoother curves than 
other methods, at the price that the curves will not exactly pass through the given 
points. The smoothness of a curve is mathematically expressed in terms of the 
continuity of its parametric representations x(t) and y(t) and their derivatives. 
B-spline curves have the property that even the second derivatives x"(t) and y"(t) 
are continuous in the points where two successive curve segments meet. In Fig. 2.18 
we can see how curves appear if their zeroth, first or second derivatives are not 
continuous in some point. Although the curve in Fig. 2.18(c) is generally considered 
smooth it does not obey the strict rules imposed by the B-spline method. 

After this informal discussion we wish to see this method in action. We shall use a 
parametric representation of curves. Any point of a curve segment between two 
successive given points P and Q will have co-ordinates x(t) and y(t), where t 
increases from 0 to 1 if the curve segment is followed from point P to point О. We 
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Fig. 2.18. (a) Oth derivative not continuous; (b) 1st derivative not continuous; (c) 2nd 
derivative not continuous 


can think of ¢ as time. If we are given the points 


Ру(хо, yo) 
Pi(xi, yi) 


Edit. nd 


then the B-spline curve segment between two successive points P; and P;,, are 
obtained by computing x(t) and y(t), where t grows from 0 to 1: 


X(t) = {(a3t + a5)t + a4)t + do 
y(t) = {(b3t + b5)t + bi}t + bo 
These equations contain the following coefficients: 
45 = (—xi-1 + 3x; — Зх; + x;,2)/6 
а= (xi-1 — 2x; + x; 41)/2 
а= (—xi-1 + xi41)/2 
do = (х: + 4X; + х;+1)/6 


(2.13) 


and b;, b2, bı, bo are derived from y;.,, у, ут, Yi+2 in a similar manner. The 
formulae above are suited for efficient computation. The value of x(t) is given 
according to Horner's rule rather than in the usual notation for a polynomial. The 
coefficients a3, a5, а}, ag are computed only once for each of the curve segments, 
which is most important, since we wish to compute x(t) and y(t) a great many times 
on a single curve segment. 

To obtain some insight into the properties of the curve in the points where two 
segments meet we can investigate the function x(t) and its first and second 
derivatives for the values t = 0 and t= 1 (y(t) having similar properties): 


x(0) = ao = (х; + 4x; + x;,1)/6 
x(1) = aş + a5 + a, + a 
Using Eqs (2.13) and simplifying we have 
x (1) = (x; + 4xi+1 + X;+2)/6 
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Fig. 2.19. Three successive points 


We observe that x(0) is not exactly equal to the x-co-ordinate x; of P;, but that it is 
also influenced by the points P; , and P;,,. In Fig. 2.19 point B is the endpoint of 
segment AB but also the startpoint of BC. From the first point of view we have 


A=P,, B-P;, С=Р+ 
хв = x(1) = (xA + 4хв + хс)/6 


where xg denotes the computed x-value at B. Viewed from segment BC, however, 
we have 


A = Pisi B=P,, C=Pi41 
xg = X(0) = (x4 +4хв+хс)/6 


We have found that both ways of calculating x give the same answer, which means 
that x(t) is continuous in B. 

Differentiating x(t) twice, we find x'(t) and x"(t). We can then substitute t = 0 and 
t = 1 in these functions as we did in x(t). It is then seen that they too are continuous 
in point B. Since y(t) and its first and second derivatives are continuous as well, it 
will now be clear that B-spline curves are very smooth. 

For any curve segment between P; and P;,, we also use the points P;_, and P;,;. 
This implies that the first curve segment will be between P, and P», and the final one 
between Р„_› and P,,_,. Thus the start and the end points of the entire curve will be 
near P, and P,,_; and not near P, and P,. The following program reads the numbers 


n 

Xo Yo 
Xi У 
Xn Yn 


from the file CURV.DAT. In the output each of the п + 1 points are plotted in the 
form of a cross. Then the B-spline curve is drawn: 


/* CURVFIT: Curve fitting using B splines */ 
#include <stdio.h> 

#define MAX 100 

#define М 30 


32 


main() 
{ float xtMAX1, yCMAX], eps=0. 04, X, Y, t. xA 
y^, yB, yC, yD, aO, al, a2, аз, bO, bl, 
int n, i, J, first 
FILE #fp; 
fpzfopen("curv.dat", "r"); 
fscanf(fp, "d", %п); 
for (1=0; i£=ni i++) fscanf(fp, "Xf Xf", x*i, 
initgr()i 
/* Mark the given points: */ 
for (1=0; it=ni i++) 
< X=xCidi Y=y[i]; 
move(X-eps, Y-eps)i draw(X+eps, Yteps)i 
move(X+eps, Y-eps); drau(X-eps, Y+eps); 
} 
first=1; 
for (1=1; itn-1i i++) 
€ xA-xLti-11; xBzxLil; xC=xCi+1]; xD=xCi+2]; 
yA=yli-1]; yB=ylili yC=yli+1]; yD=yli+2]; 


a3=(-xA+3#(xB-xC)+xD)/6. Oi 

a2=(xA-2#xB+xC)/2. Oi 

al=(xC-xA)/2. Oi 

a0=(xA+4#xB+xC)/6. Oi 

for (J=0; J<=Ni g++) 

< t=(float) y/(float)Ni 
X-2CCGagstea2)5tal)sttaO; 
Y=((bS#t+b2)#t+b1)#t+bO0; 
if (first) < first-O; move(X, 

y 


bi=(yC-yA)/2. Oi 


} 


endgr(); fclose(fp)i 


› 
The following data 


.85 


сл 


л CA 


N 
CA 
к= ке ка ка Ба ка ка ка ка к 
ьо ә н NN 


SS. 46 
in file СОКУ. DAT will produce the output of Fig. 2.20. 


Yi} else draw(X, 
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xB, 
b2, 


xC, 
b3; 


xD, 


у+1); 


b3=(-yA+3#(yB-yC)+yD)/6. Oi 
b2-(yA-2*yB*yC)/2.0; 


bO=(yA+4#yB+yC)/6. О; 


т? 
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Fig. 2.20. Output of program CURVFIT 


EXERCISES 
2.1 Write a program to draw a set of N ellipses whose parametric equations are: 
x = xo + (iR/N) cos ф 
y= уо + ((N — i)R/N)sin p 


Choose fixed values xo, yo, R, N; for example, xo = 4, у= 3.5, К =3, N = 40. 


Let i range from 1 to N — 1. For each ivalue, let ф successively have the values 
0°, 6°, 12°,... , 360°. 


Fig. 2.21. Tree 
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2.2 


2.3 


2.4 


2.5 
2.6 


2.7 
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Write a program to draw a set of 30 triangles. The vertices of the first triangle 
are (1, 1), (6, 1), (3, 5). Each next triangle is obtained by rotating the previous 
one about point (4, 2) through an angle of 3°. 

Implement a recursive line-clipping algorithm based on bisection. Suppose a 
window and a line segment РО are given. If P and О are both inside the 
window, РО can be drawn. This is also the case if one of these points is inside 
the window and the other is on a window edge. There are also some cases 
where we can easily decide that nothing is to be drawn. Invent these yourself. 
In all other cases, find point M in the middle of PO and apply the same 
procedure recursively to PM and МО. Choose a tolerance with respect to the 
window edges, and investigate the influence of this tolerance on the number of 
bisections. 

Develop and implement an algorithm for clipping a line against a window 
which is a triangle. 

Modify program STARS in Section 2.7 such that stars will not overlap. 

Use random numbers to generate a sequence of, say, 30 points, and apply 
curve fitting to draw a curve (approximately) through these points. 

Write a program which produces a more realistic tree than Fig. 2.17. An 
example is given in Fig. 2.21. 
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Geometric tools for three-dimensional 
algorithms 


3.1 VECTORS 


A certain amount of mathematical knowledge is essential to understand and write 
programs for three-dimensional graphics. This book is not intended to be a textbook 
on mathematics, and the reader is assumed to be already familiar with the 
mathematical topics of this chapter, especially with vectors and determinants. 

A vector is a directed line segment, characterized by its length and its direction 
only. Figure 3.1 shows two representations of the same vector а= PQ = b = RS. 
Thus a vector is not altered by a translation. In Fig. 3.2 the initial point of b is the 
terminal point of a. Then the sum of a and b is defined as the vector c drawn from 
the initial point of a to the terminal point of b, and we write 


c=a+b 
The length of a vector a, denoted by |a|, is the distance between its initial and its 
terminal points. A vector with length zero is the zero vector, written 0. The notation 
—a is used for the vector that has length |a| and whose direction is opposite to a. For 
any vector a and real number c, the vector ca has length |c] |a|. If a=0 or c —0, 
then ca = 0, otherwise ca has the direction of a if c > 0 and the opposite direction if 
€ « 0. For any vectors u, v, w and real numbers c, k we have: 


u+v=v+u 

(u+v)+w=u+(v+w) 

u+0=u 

u+(—u)=0 

c(u + v) = си+су 

(c+k)u=cu+ku 

c(ku) = (ck)u 

lu=u 

Ou = 0 

Figure 3.3 shows three unit vectors i, j, k. They are mutually perpendicular and 

have length 1. Their directions are the positive directions of the co-ordinate axes. 
We say that i, j, k form a triple of orthogonal unit vectors. The co-ordinate system is 
right-handed, which means that if a rotation of i into the direction of j through 90° 


corresponds to turning a right-handed screw, then k has the direction in which the 
screw advances. 
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[д] 


о; 


Fig. 3.1. Equal vectors 


O; 
oF 


a 
Fig. 3.2. Vector addition 
We often choose the origin O of the co-ordinate system as the initial point of all 
vectors. Any vector v can be written as a linear combination of the unit vectors i, j, 
k: 
v=xi+yj+zk 
The real numbers x, y, z are the co-ordinates of the terminal point P of vector 
v= OP. We also write this vector v as either a row or a column: 
X 
у= [х у 2] or у= |у 
2 


The numbers х, у, z are sometimes called the elements of vector v. 


y 


Fig. 3.3. Right-handed co-ordinate system 
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3.2 INNER PRODUCT 


The inner product, dot product or scalar product of two vectors a and b is written 
а · b and is defined as 
a-b=|a||b|cosy if а+0 and b+0 


. (3.1) 
a-b=0 if a=0 or b=0 
where y is the angle between a and b. 
Applying this to the unit vectors i, j, k we find 
i-i=j-j=k-k=1 
ы (3.2) 


i-j=j-i=j-k=k-j=k-i=i-k=0 
Setting b =a in Eqs (3.1) we have a: a=|al*, so 
|а| = Vla - al 
Important properties of inner products are 
c(ku * v) = ck(u * v) 
(си + Ку): = сиу + kv: w 
u:v=v-u 
u-u=0 опу u=0 
The inner product of two vectors u = [u, из из] and v = [v, v; v3] can be computed 
E U’ V= иу + u5U5 t изо; 
This is proved by writing the right-hand side of 
u‘ v = (иі + и) + uk) * (vii + vaj + vsk) 


as the sum of nine inner products and then applying Eqs (3.2). 


3.3 DETERMINANTS 


Before proceeding with vector products we shall pay some attention to 
determinants. 
To solve the following system of two linear equations 


ese 63 
we multiply the first equation by b>, the second by —b, and add, finding 
(a,b, — a5bi)x = bsc, — Бус» 
Then we multiply the first equation by —a5, the second by a, and add again, finding 
(a,b, = a5b,)y = ас — anc, 
If a,b, = a,b, is not zero, we can divide and find 


_ b2€ı > bic; _ 44,€C5 — 45€, 


= ; = 3.4 
a,b, — ab, ab; = ab, (3:9) 
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The expression in the denominators is written in the form 


a, bi 


аз b, 
and then called a determinant of second order. Thus 


a, b, 
= a1b, а,б! 


а Б, 
With determinants the solution of Eqs (3.3) can be written 


D, D; 


=—, =—, (D+0 
imr pes Qe 
where 
D- a, bj Bos C1 bi D, = a, €i 
аз b, €; b; аз C2 


Notice that D, is obtained by replacing the ith column of D with the right-hand side 
of Eqs (3.3) (1= 1 or 2). This method of solving a system of linear equations is 
called Cramer's rule. It is not restricted to systems of two equations (although it 
would be very expensive in terms of computer time to apply the method to large 
systems). We define determinants of third order by the equation 


b 
i di b, c; b, с, b, c 
D- d2 b, C2 а, — 05 + 43 
b3 сз b3 сз b, c; 
аз b3 c3 
and determinants of fourth order by 
ау b, €1 d, 
b d 
Dos a2 2 C2 2 
аз b3 сз d; 
ал b, C4 d, 
b; C2 d, b, Ci d, b, Cy d, b, Ci а, 
= а, b; C3 d; — d2 b; C3 d; + a3 b; C2 d; = 44 Б, C5 d; 
b, cy di by Ca di; bj €4 dy Б; сз d 
апа $о оп. 


Determinants have many interesting properties, some of which are listed below. 


(1) The value of a determinant is not altered if its rows are written as columns 
in the same order. For example: 

a, bi 

a, b, 


a, а 


b, b 


(2) If any two rows (or two columns) are interchanged, the value of the 
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determinant is multiplied by —1. For example: 


а, b, с, ay by с, 
а, b; о|=- |a; b3 сз 
аз b; C3 a; b; C5 


(3) If any row (or column) is multiplied by a factor, the value of the 
determinant is multiplied by this factor. For example: 


a, b, 
аз b, 


ca, cb, 
a, b, 


(4) If a row (or a column) is altered by adding any constant multiple of any 
other row (or column) to it, the value of the determinant remains 
unaltered. For example: 


а, b, с, а, bi с, 
а b; C2 =|а; b; c; 
аз + ka, b;+kb, сз + Кс, аз Ёз сз 


(5) If a row (or column) is a linear combination of some other rows (ог columns), 
the value of the determinant is zero. For example: 


а, b, Ci 
а b; C2 es () 
3a, = 2a> 3b, == 2b, 3c, == 2с› 


There are many useful applications of determinants. Determinant equations 
expressing geometrical properties are elegant and easy to remember. For example, 
the equation of the line in R, through the two points P,(x;, yı), P2(x2, у) can be 
written 


x y 1 
Х| y 1 = 0 (3:9) 
X2 у 1 


This can be understood by observing first, that Eq. (3.5) is a special notation for a 
linear equation in x and y, and consequently represents a straight line in R;, and 
second, that the co-ordinates of both P, and P; satisfy this equation, for if we write 
them in the first row we have two identical rows. Similarly, the plane in К; through 
the three points P(X, у, 21), Po(x2, у, 22), P3(x3, ys, z3) has the equation 


X y g 1 
xy yy zn 1 
о y 25 1 

1 


Хз Уз 23 


=0 
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3.4 VECTOR PRODUCT 
The vector product or cross-product of two vectors a and b is written 
axb 


and is a vector v with the following properties. If a — cb for some scalar c, then 
v — 0. Otherwise the length of v is 


Iv| = lal |b| sin y 


where у is the angle between a and b, and the direction of v is perpendicular to both 
a and b and is such that a, b, v, in that order, form a right-handed triple. The latter 
means that if a is rotated through an angle y < 180? into the direction of b, v has the 
direction of the advancement of a right-handed screw if turned in the same way. The 
following properties of vector products follow from this definition: 


(ka) xb  —k(aX b) for any real number k 
aX(b+c)=axXb+axXe 
aXb =—bXa 


In general, aX (bX c) £ (aX b) X c. 
Using a right-handed orthogonal co-ordinate system as in Section 3.1, with unit 
vectors i, j, К, we have 


ixi=jxj=kxk=0 
ixj=k, jx k=i, kXi=j 
}хї=—К, kXj- -i, iXk- -j 
Using these vector products in the expansion of 
a X b = (аі + a,j + ask) X (bj + bj + bk) 
we obtain 
a X b = (ab; — a3b2)i+ (ab, — a,b3)j + (abo — ab ,)k 


which can be written 


а а; |, аз а, а az 
ахь = + + К 
b; b; b, b, b, b; 
We rewrite this in a form that is very easy to remember: 
i J È 
aXb-j|a, а, а; 
b, b, bi 


This is a mnemonic aid rather than a true determinant, since the element of the first 
row are vectors instead of numbers. If a and b are adjacent sides of a parallelogram, 
as in Fig. 3.4, the area of this parallelogram is the length of vector a X b. This 
follows from |a X b| = |а| |b| sin y. 

In Fig. 3.5 the vectors a and b lie in the plane through the x- and y-axes. We 
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E 


= 


Fig. 3.4. Parallelogram with area |a X b| 


imagine a z-axis whose positive direction is out of paper towards the reader, so the 
co-ordinate system is right-handed. Then in three dimensions we have 


а = [a, d2 0], b = [b, b; 0] 


i j 
а, а 
aXb- |а, a, 0|-|, ML 
b, b, 0 dian: 


Thus the product vector aX b has the same direction as k if and only if the 
determinant 


is positive. This implies that the rotation of a into b through an angle less than 180? 
is counter-clockwise if and only if D >0. We shall use this principle to determine 
whether the vertices A, B, C of a triangle, in that order, are traversed 
counter-clockwise. 


Fig. 3.5. Vector product k =a X b 
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ct 


Fig. 3.6. A, B, C counter-clockwise 
In Fig. 3.6 we have 
u-[u, и] = AB, v=[v,; v,]=AC 


XA Ya 1 XA YA 1 
D-|xs ys 1|=|Ххв—Хд ув—уд 0 
1 


Xe Ye Xc—XA ус—уд 0 
_ |XB^XA Ув-УА|= | и U2 
Xc—XA yc-YA!—|U U2 


The vertices A, B, C, in that order, are traversed counter-clockwise if and only if u 
is rotated into the direction of v through an angle y < 180° counter-clockwise. This 
means that we can derive the rotation sense of A, B, C from the determinant 


ta Ja 1 
D= | XB ув 1 
xe yc 1 


in the following way: 


If D >0, the points A, B, C, in that order, are traversed counter-clockwise. 
If D <0, the points A, B, C, in that order, are traversed clockwise. 
If D =0, the points A, B, C lie on the same straight line. 


3.5 DECOMPOSING POLYGONS INTO TRIANGLES 


In Chapters 4 and 5 we shall produce pictures of three-dimensional objects whose 
boundary surfaces will be polygons. This is no serious restriction since curved 
surfaces can be approximated by a great many polygons in the same way as a curve 
is approximated by a sequence of line segments. Dealing with arbitrary polygons can 
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Fig. 3.7. Two polygons partly hiding each other 


lead to quite complex situations, especially if visible and hidden line segments are 
distinguished. Figure 3.7 shows an example of such a situation. 

If all interior angles of a polygon are less than 180? the polygon is called convex. 
In Fig. 3.8(b) the interior angle at vertex P is greater than 180°. We shall call such a 
vertex non-convex. All other vertices in Figs 3.8(a) and 3.8(b) are convex. If a 
polygon has at least one non-convex vertex, the polygon itself is said to be 
non-convex. 

If A and B are two points on a convex polygon the entire line segment AB 
belongs to the polygon. For a non-convex polygon this may not be the case. 
Non-convexity of polygons is a source of complexity and so is their variable number 
of vertices. For these reasons we pay special attention to triangles. These obviously 
have a fixed number of vertices and they are always convex. They are also 
interesting in connection with arbitrary polygons because any polygon can be 
divided into a finite number of triangles. This will be the subject of this section. 

Division of a convex polygon into triangles is extremely simple, as Fig. 3.9(a) 
shows. If the vertices are successively numbered Po, P,,...,P,_,, then drawing the 
diagonals PoP», PoP3, . . . , PoP„—2 is all that is needed. In a non-convex polygon such 
as in Fig. 3.9(b), however, this simple method will not work, since some of the 
diagonals РР», PoP3,...,PoP,-2 may not completely lie inside the polygon. We 
shall now develop a program which reads the co-ordinates of a polygon's vertices 


(a) (b) 
Fig. 3.8. (a) Convex polygon; (b) non-convex polygon 
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«О 


oU 


(a) (b) 
Fig. 3.9. (a) Diagonals inside polygon; (b) diagonal P,P; not usable 


and perform a division into triangles.' We shall require the vertices to be specified in 
a counter-clockwise order. For example, in Fig. 3.9(b) the sequence P,PsPsPoP,P2P3 
will do, but the sequence P4P:P,P;P;P,P, will not. For a polygon with n vertices the 
number п is given first, followed by the п co-ordinate pairs of successive vertices in 
counter-clockwise order. The output will be a drawing of the polygon in which 
diagonals are drawn, completely dividing the polygon into triangles. Before drawing 
the diagonal we have to ensure that the entire diagonal lies inside the polygon. 
Suppose that P; ,, P; P;,, are three successive vertices, where we define 
P ,—P, , and P, = Po, to allow the cases i = 0 and i = n — 1. Remember that the 
vertices are given in counter-clockwise order. Then P; is a convex vertex if and only 
if the three vertices P; ;, P; P;,,, in that order, are also traversed counter- 
clockwise. As a counter-example, consider Fig. 3.10, where triple P,P2P; is 
clockwise, P, is non-convex and diagonal P,P, lies outside the polygon. Thus 
diagonal P; ,P;,, can only be a candidate for our purposes if P; (x; 1, Yi-1) 
P,(x;, у), Р+1(Х+1› yiii), in that order, are traversed counter-clockwise, that is, if 


Xii Ver d 
D=| x y 1|>0 
Ху Jii 1 


Po P 


Fig. 3.10. Diagonal P, P, outside polygon 


! We shall use a method which works properly for a large class of polygons. However, there 
are polygons for which the method will fail, see Exercise 3-3 at the end of this chapter. 
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Fig. 3.11. Diagonal P,P, partly outside polygon 


This condition is necessary, but, unfortunately, it is not sufficient, as is shown by 
Fig. 3.11. Неге Po, P;, P2, in that order, are traversed counter-clockwise yet PoP» 
cannot be used to divide the polygon into triangles. This situation can be avoided by 
taking the lengths of the diagonals into account. We shall choose the shortest 
diagonal P;_; P;,; that has a convex vertex P; between P; , and P;,,. This diagonal 
is used to cut off triangle P;. , P; P;,,. Then the remaining polygon 


Ре, Pis و‎ Pia Piri +» Р 
is treated in the same way, and so on. Technically this is realized by introducing an 
integer array vo,..., v,.,, containing the vertex numbers of the remaining 
polygon. Initially we set m =n and v; =i (i= 0, 1, ..., n — 1). Every time a triangle 


is cut off, m is decremented. If a program like this is meant for practical use many 
tests on the validity of the input data are desirable. A program that adequately 
rejects any set of invalid input data is said to be robust. In our situation, it makes 
sense to test whether the given input sequence 


п Xo Yo Xi у... Xn-1 Yi 
describe a polygon at all. This is not the case with the sequence 
41122 21 12 


since traversing the points in the given order will yield Fig. 3.12, which we do not 
accept as a polygon. 
Other tests that are appropriate concern 


The maximum value of n, e.g. n < 500; 
The minimum and maximum value of the co-ordinates; 
The required counter-clockwise orientation. 


Despite their importance, most of these tests have been omitted here and left as an 
exercise for the reader. On the other hand, the program will contain a special 
feature, which, strictly speaking, could also have been omitted, namely the 
representation of diagonals by dashed lines instead of by unbroken lines. We shall 
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Fig. 3.12. Result of invalid sequence 


require all dashes of a dashed line to have the same length. A dashed line must not 
start or end with a gap, but there should be a full dash at the beginning and at the 
end of the dashed line, as shown for segment PO in Fig. 3.13. 


P Q 


Fig. 3.13. A dashed line 


The reader is encouraged to solve this problem and to compare his solution with 
the function dash in the following program. 


/* POLY, TRIA: Dividing a polygon into triangles */ 
define NMAX 500 

define BIG 1.0e30 

#include <math. h> 

int n, vENMAX]; float хСММАХЈ, yCNMAX Ji 


main() 
€ int i, h, ys m ls imini 
double diag, min diag; 
printf("4sNnZsNn", 
"Give n, followed by n coordinate pairs (х, у) of the vertices, ", 
"in counter-clockuise order"); 
scanf("%d", п); if (n>=NMAX) error("n too large"); 
for (i20; іп; i++) < scanf("Xf Xf", £xLil, £yLil)i vlil=i; } 
initgr(); draw polygon(); т=п; 
while (m>3) 
< min diag-BIG; 
for (1=0; 1<т; i++) 
€ h= (i==0 ? m-1 : 1-1); у= (i==m-1 ? О: 1+1); 
if (counter clock(h, i, J, &diag) && diag<min_diag) 
€ min diag-diag; imin=i; 
+ 
+ 
i=imini h= (і==0 ? m-1 : 1-1); y= (i--m-1 ? О: i*1)i 
if (min diag--BIG) error("urong sense of rotation"); 
dash(xCvCh1], y£vCh1], xtvELj13, ytvEj12); 
m--i 
for (l-i; lm; l**) v[l]=v[l+1]; 


DECOMPOSING POLYGONS INTO TRIANGLES 


endgr(); 
3 


error(str) char stri 
< printf("ZsNn", str); exit(1); 
3 


int counter clock(h, i, J, pdist) int h; і, Ji double #pdist; 


{ double xh-xtvEh11, xi-xtvEtill, xy=xCvl 1], 
gh-ytvEh11J, yi=yCvCill, yj-yLvEj11, 
x hi, gy hi, x hj; y hy, Determi 
x hi-xi-xh; y һі=уі-уһ; x hj-zxj-xh; y hJ-yJ-yh; 
#pdist = x hj + x hj + y hj * y hyi 
Determ = x hi * y hj — x hj * y hii 
if (Determ == О) 


error("Three successive vertices on the same straight line"); 


return Determ»Oi; 
} 


draw _polygon() 
< int ii 

move(xtn-11, уСп-11); 

for (1=0; i£n; i++) drau(xLil, ytil); 
} 


dash(x1, yi, x2, y2) float xi, yi, x2, y2i 

€ int i. ki 
float xdif=x2-x1, ydif=ye-yl, pitchO=0.3, dx, dui 
k = 2 * ceil(sqrt(xdif*xdif*ydifsydif)/pitchO) + 1; 
dx=xdif/ki dy=ydif/ki 
for (1=0; i£k; i*-2) 
€ move(xiti#dx, уі+ізау); drau(xi*(iti)sdx, y1+(Ci+1)#dy); 
> 


The following input sequence 


Fig. 3.14. Output of program POLY-.TRIA 


47 


48 GEOMETRIC TOOLS FOR THREE-DIMENSIONAL ALGORITHMS 


3.6 HOMOGENEOUS CO-ORDINATES 


This section deals with a mathematically interesting topic which is related to 
perspective and therefore often included in books on three-dimensional computer 
graphics. We shall conform to this usage in order to obtain a better understanding of 
notations such as [x у 1]апа[х y z 1]. This section deals with geometry, and 
though not particularly difficult, it is rather long and theoretical. Readers who are 
interested primarily in practical aspects of graphics programming might skip to 
Section 3.7 without serious consequences for the study of the rest of the book. They 
can return to this point at a later stage. 

In Section 2.3 we used the notation [x y 1] to denote a matrix of only one row, 
sometimes called a row vector. This notation could also have been introduced as a 
special case of [x y w] where x, y, w are called homogeneous co-ordinates. We 
use three homogeneous co-ordinates to denote a point in two-dimensional space 
(2-space, for short). In projective geometry homogeneous co-ordinates had been 
used long before computer graphics became popular. In Chapter 4 we shall deal 
with perspective transformations. They are in fact central projections, treated 
thoroughly in books on projective geometry. We shall briefly discuss only a few of 
the topics of this fascinating but rather difficult branch of mathematics, avoiding 
formal definitions and rigorously proven theorems. 

In Fig. 3.15 we have an x-axis and a w-axis, so a point is given by its co-ordinate 
pair (x, w). Any point P(x, w) not on the x-axis has its central projection P'(X, 1), 
being the intersection of line OP and line / with equation w — 1. The origin O is the 
centre of projection. Line segment PO can be considered as a ray of light from an 
object P to the eye O. Introducing the points O(0, w) and Q'(0, 1), we have two 
similar triangles ОРО and OP’Q’, so that 


.X PQ PO x 
`1 OQ' OQ w 


0 x/wW X x 


Fig. 3.15. Two-dimensional central projection 
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All points (x, w) with the property x = wX lie on line OP and have the same 
projection P'. If we are interested only in projections on line /, not the actual values 
x and w but only their ratio matters. It is natural to use only one co-ordinate X 
instead of a co-ordinate pair (X, 1) if we deal only with points on line /. Should we 
insist on using a co-ordinate pair, any number pair (x, w) satisfying x/w — X would 
do if we so agreed. In geometric terms the co-ordinate pair (wX, w) of any point P, 
distinct from O, on line OP' could serve as a notation for point P'. This is what is 
done when homogeneous co-ordinates are used. In general, any point 
(Xi, X5, ..., X,) in n-space is written as a point (WX, wX5,..., wX,, w) in 
(n + 1)-space, where w is any non-zero real number. The latter n + 1 numbers are 
said to be homogeneous co-ordinates of the original point in n-space. The reader 
will probably be familiar with the notion of projection, being a many-to-one 
mapping from (n + 1)-space to n-space. Homogeneous co-ordinates arise from the 
reverse process, a one-to-many mapping from n-space to (n + 1)-space. Point (5, 7) 
in 2-space, for example, can be written in homogeneous co-ordinates as (15, 21, 3), 
or as (500, 700, 100), etc. Though mainly interested in 2-space, we can consider 
these triples as points in 3-space. Obviously, (X, Y) is the usual 2-space notation for 
point (X, Y, 1) in 3-space. This point, P', is the central projection of any point 
P(x, y, w), if x/w — X and y/w — Y. Again the origin O is the centre of projection; 
all point are now projected to the plane w = 1. 
To explain the term homogeneous let us consider the equation 


aX +bY+c=0 (3.6) 
representing a line in 2-space. Replacing X and Y with x/w and y/w, we have 


a(x/w) + b(y/w)+c=0 
iin ax + by +cw =0 (3.7) 
It is customary to call Eq. (3.7) homogeneous, because of the identical structure of 
the terms ax, by, cw. It is therefore reasonable to call x, y, w homogeneous 
co-ordinates of point (X, Y). If the 2-space is again the plane w=1 in an 
xyw-co-ordinate system, then Eq. (3.7) is the plane through the origin O and the 
given line. 

If (x, y, w) were merely used as another notation for (x/w, y/w) it would be 
necessary to stipulate that w be non-zero. However, in this way homogeneous 
co-ordinates would hardly have any advantages over conventional ones, and the 
reader will probably wonder whether they have any at all. To do them justice, we 
have to admit that w — 0, as we shall see presently. Homogeneous co-ordinates lead 
to elegant and general formulations of geometrical properties. For example, let us 
consider a system of two linear equations, each representing a line in 2-space: 


ане ае 


aX + bY +c = 0 (3.8) 


If е two lines are parallel, they have no point of intersection and there is no 
number pair (X, Y) satisfying Eqs (3.8). Thus to find a common point of the two 
lines we have to use a rule with an exception, which is not particularly elegant. 
Replacing X and Y with the homogeneous co-ordinates x, y, w will improve the 
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situation: 


ears dicata (3.9) 
ах + boy + caw = 0 


Equations (3.9) can also be interpreted as planes through O; they have at least the 
trivial solution x = y = w —0. For a geometrical interpretation we choose concrete 
values for the coefficients. Let us, for example, replace Eqs (3.8) with 
ind -t3Y- 6-0 
4X + 6Y – 24=0 


representing the parallel lines of Fig. 3.16. We then have to replace Eqs (3.9) with 


[os би = 0 (3.10а) 
4x + бу – 24w = 0 (3.106) 
This system is equivalent to 
(2 + 3y = 0 
w = 0 


so the solution consists of all triples (3k, —2k, 0), where К is any real number. In 
3-space these points constitute the line through O and (3, —2, 0), this line being the 
intersection of the planes given by Eqs (3.10a, b). Returning to the 2-space of the 
plane w=1, we remember that every point (X, Y) is associated with the line 
(wX, wY, w) in 3-space. For non-zero w-values this association is almost trivial. An 
important motivation for the use of homogeneous co-ordinates will now become 


Fig. 3.16. Parallel lines 
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clear. To every line in 2-space we add an object called infinite point. This infinite 
point cannot be written in conventional co-ordinates, but it can in homogeneous 
co-ordinates. For example, the infinite point of the line of Eq. (3.10a) is written as 
(3, —2, 0) or as any other triple (3k, —2k, 0) with non-zero К. Since these triples are 
solutions to both Eqs (3.10a) and (3.10b) this infinite point is said to be the 
intersection of the two parallel lines in Fig. 3.16. Note that this infinite point of 
intersection is considered to be a point in 2-space. As we have seen, every point in 
2-space is associated with a line in 3-space, so we might wonder what line the infinite 
point (3, —2, 0) is associated with. Since this has to be the line through the given 
point and O(0, 0, 0), the line we looked for is the one through O and (3, —2, 0), 
lying in the plane w =0. 

It is reasonable to call (3, —2, 0) an infinite point, since we can regard it as the 
limit of (3, 2, и) where w goes to 0, and the latter triple in homogeneous 
co-ordinates is equivalent to (3/w, —2/w, 1) which lies far away for small values of 
и. The notion of infinite points enables us to say that any two distinct lines in 
2-space meet in one point. In the same way we say in projective geometry that any 
two distinct planes in 3-space have a line of intersection. If the planes are parallel, 
all points on this line of intersection are written (x, y, z, 0) in homogeneous 
co-ordinates. We shall not discuss this in greater detail, but return to 2-space and 
show other new possibilities offered by homogeneous co-ordinates. 

With non-homogeneous co-ordinates a linear transformation in 2-space can be 
written 


[Х' Y']=[X Y]A 

a, a2 

5 

bı b, 
Since [1 0]A = [а, a] and [0 1]A=[b, bz], the rows of matrix А are the 

images of [1 0] and [0 1], respectively. 

No matter how A is defined, the origin O is self-corresponding, that is, 
[0 0]A = [0 0], so a translation cannot be expressed in this way. In homogeneous 


co-ordinates, however, a point in 2-space is given by a triple (x, y, w) and a 
transformation is written 


k y" w']=[x у wJA 


dı а аз 
А = b, b, b; 
Су С C3 


We have 
[1 0 OJA=[a, а, aj] 


Point [1 0 0] is the infinite point of the x-axis. Thus the first row [a a, as] of 
matrix A is the image of the infinite point of the x-axis. Similarly, the second row 
[bı b2 bı] is the image of the infinite point of the y-axis. Since 


[0 0 1]А=[с, с, ез] 


we see that the third row [су с: сз] is the image of the origin [0 0 1]. This 
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B (by, b», 1) 
y (01,0) Í Fá 
| ^ 
/ 
| ў 
| x 
bi 
(0,1,1) 4“ 
0 (0, 0, 1) P(1,0,1) x Clcyc2,1) P' A(a,,a5, 1) 


Fig. 3.17. Quadrant mapped into a triangle 


means that homogeneous co-ordinates enable us to express any translation by a 
matrix multiplication. In fact this is nothing new, since we already used it in Section 
2.3. However, translation is not the only new possibility offered by this type of 
matrix multiplication; we can also map parallel lines to intersecting lines. We shall 
show this, mapping a complete rectangular quadrant to a triangle. 

In Fig. 3.17 we have an arbitrary triangle whose vertices are in rectangular 
co-ordinates A(a,, a2), B(bi, b2), C(ci, c2). We have added a third co-ordinate 1 as 
a formal means of obtaining homogeneous co-ordinates. Exceptions are the infinite 
points (1,0, 0) and (0, 1,0) of the co-ordinates axes, whose third homogeneous 
co-ordinates are zero. We shall now construct a matrix M, such that the matrix 
multiplication 


kK у' 2']=[х y z]M 
maps point O to C, point (1,0,0) to A and point (0, 1,0) to B. For any non-zero 
values a, В, y the points A, B, C will be the required image points if 
aa, са, а 
M-| Bb, Bb, В (3.8) 
Ya. Үс Y 
This is easily verified, since we have, for example, 
[1 0 ОМ=[аа, aa, а] 


whose right-hand side is just another notation for point A in Fig. 3.17. At first sight 
it seems that the constants а, B, y in Eq. (3.8) are not relevant and can be set to 1. 
We need them, however, to be able to prescribe the image U' of the so-called unit 
point U. Point U' may be chosen anywhere inside triangle ABC. Since U(1, 1, 1) 
maps to U'(u,, u2, 1), we have 


[1 1 1]M=[u, u 1] 
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which is short-hand for the following system of three linear equations in а, В, y: 
aa +b,B+c;y=u: 
аа + b; + су = us 
a+ B+ y=1 


This system has a unique solution, which follows from the fact that (aj, a5, 1), 
(bi, b2, 1), (ci, сг, 1) denote vertices of a triangle. Solving this system for a, В, y 
and substituting the results into Eq. (3.8) gives us the desired matrix M. 

Straight lines map to straight lines, but parallelism is not preserved. For example, 
the vertical line through point U in Fig. 3.17 maps to the line through B and U’. 
This provides a geometrical means to find the projection P' of point P. Having 
computed matrix M we can find P' analytically as the product 


[1 0 цм 


Note that every infinite point maps to a point on AB, and that the infinite points of 
parallel lines map to a single point on AB. It makes sense to view the entire 
quadrant as a flat landscape of which triangle ABC is a picture and AB is the 
horizon. 

This concludes our discussion on homogeneous co-ordinates. For many more 
interesting properties the reader is referred to a textbook on projective geometry 
such as Hopkins and Hails (1953). 


3.7 THREE-DIMENSIONAL TRANSLATIONS AND ROTATIONS 


If every point P(x, у, z) is mapped to a point P'(x', y', 2'), according to 


x'=x+a, 
у'=у+а, 
z'=z+a; 


where а), a2, аз are constants, we have a translation in three-dimensional space. 
This translation can also be written 


[ху z 1]=[x y z ЦТ 


1.000 
0 1 0 0 

T= À 
0 0 1 0 em 


а, @ аз 1 


Readers who have studied Section 3.6 will know that the first, the second and the 
third rows of matrix T denote (images of) the infinite points of the co-ordinate axes, 
and that the fourth row is the image of [0 0 0 1]. The latter means that in 
homogeneous co-ordinates [a, a; a, 1] is the image of the origin 0. 

Rotations about the co-ordinate axes can be expressed with matrices without 
using homogeneous co-ordinates. For the sake of brevity we shall do so now and 
switch to homogeneous co-ordinates when we actually need them. We use a 
right-handed co-ordinate system, and we call a rotation about an axis positive if it 
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z 


dD 
"i 


x 


Fig. 3.18. Positive rotations about axes 


corresponds to the positive direction of that axis in the sense of a right-handed 
screw. This is shown in Fig. 3.18. 

We shall rotate about the z-axis through an angle œ, and abbreviate cos а = c, 
sin а = s. The rotation is then written: 


, , 


[I y z]s[x y AR 


e $ 0 
R,=| -s с 0 
0 0 1 


which follows from Section 2.3, Eq. (2.4). 

This matrix R; can be used to derive the matrices R, and R, for rotations about 
the other two axes in a formal way, that is, without using a picture. This is done by 
the cyclic permutation obtained by replacing each of the letters x, y, z with its 
successor, the successor of z being x. 

We convert К, into К, by cyclic advancing each row one position, followed by a 
similar operation on the columns: 


G $ 0 0 01 
R,=|-s c 0 P é $ 0 
0 0 1 WX 

1 0 0 
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In the same way R, is converted into its ‘successor’ Ry: 


1 0 0 0 == с 
R,=|0 c s =: 0 0 
0 -s c 0 c s 


Summarizing, we have the following matrices: 


1 0 0 
R,-|0 cosa sin« (3.10) 
0 -—sin« cosa 


cos æ 0 -sin« 
R,= 0 1 0 (3.11) 
snæ 0 cosa 


cosæ sina 0 
К,= | –5іп а cosa 0 (3.12) 
0 0 1 
Matrix R, is used in the following way: 


[x' y z'l[x у 2], 
for a rotation about the x-axis through an angle а, and the matrices R, and R, are 
used similarly. 

As explained in Section 2.1, equations for a transformation can also be 
interpreted as a change of co-ordinates. Moving a point a certain distance to the 
right requires the same equations as moving the co-ordinate system the same 
distance to the left. In practice it is more convenient to move the co-ordinate system 
in the same sense as points would have been moved, but this requires inverted 
matrices. Fortunately, the inverses of T, Rx, Ry, R, (Eqs (3.9)-(3.12)) can be 
written down immediately: 


1 0 0 0 
0 1 0 0 
T t= 3.13 
0 0 1 0 ( ) 
= — 4; — 43 1 
1 0 0 
R,'=|0 cosa —sin « (3.14) 


0 sina соѕа 
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cosa 0 sin« 
Ry = 0 1 0 (3.15) 
—sinæ 0 cosa 


cosa -—sin« 0 
R;'-|sina cosa 0 (3.16) 
0 0 1 


We are now in a position to find the matrix К for a rotation about any line passing 
through the origin O. To define the sense of rotation we prefer to say that we rotate 
about a vector v whose initial point is O. Then a positive rotation corresponds to the 
vector direction according to a right-handed screw. As before, we shall rotate 
through an angle a. 

If the endpoint of v is given in rectangular co-ordinates we first compute its 
spherical co-ordinates р, Ө, q (see Fig. 3.19): 


p = | = V(vi + v3 + v3) 
If p =0, we set 0 = ф = 0. Otherwise 


arctan (v;/v;) if v,>0 

e- л + arctan (05/0) if v,«0 
z/2 if v,-0 and у, 20 
32/2 if v,;=0 and v,<0 


Ф = arccos (vs/p) 


Fig. 3.19. Spherical co-ordinates 
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(The reader will probably be more familiar with the reverse computation: 
v, =p sin ф cos 6, v= p sin ф sin Ө, u; = p cos ф) 


Our strategy will now be to change the co-ordinate system such that v, the axis of 
rotation, will lie on the new positive z-axis. We begin with a rotation of the x- and 
y-axes about the z-axis through an angle Ө. According to Eq. (3.16) we can write 
this as 


' 


kw y' :']=[х у z]RZ 


cos0 -sin 0 
R;'-|sinÓ cos0 0 
0 0 1 


The x'-axis has the positive direction of vector (v, v; 0). Then we rotate the x'- and 
the z'-axes about the y'-axis through an angle Ф, yielding a z'"-axis in the direction of 
vector v (see Fig. 3.19). 

Referring to Eq. (3.15) we write this as 


[x" y" z"]  [x' x' z']R;" 


cosp 0 sing 
R,'- 0 1 0 
—sing 0 соѕф 


The actual rotation about v through an angle a can now be performed as а 
rotation about the z’-axis. Referring to Eq. (3.12) we have 


m m 


[x y z"] = [х” у" z"|R, 


cosa Sina 0 
R,-|-sin« cosa 0 
0 0 1 


So far, we have achieved the following: 
x" y" a] = [х y 2]К; 'К;'К, 
Unfortunately, the co-ordinates х”, у”, 2" refer to the latest co-ordinate system, 
whereas we want them expressed as the original co-ordinates. We shall denote these 


original co-ordinates of the rotated point by x*, y*, z*. Reverting to the original 
co-ordinates is accomplished by applying the inverses of matrices R;' and Ry 


т m m 


(being R, and R,) in reverse order to x”, y”, z^: 


[x* y* z*] — [x" y" z"|R,R, 
This means that the complete rotation about vector v through the angle a is 


computed as follows: 
[х* y* «*]=[х y z]Rz'R;'R,R,R, 
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cos@ -—sinO 0 
R;'-|sinÓ cos@ 0 


0 0 1 


where 


cog 0 sing 
0 1 0 


لا 


—sing 0 со$ф 


cosa sina 0 
R= E а cosa ^ 
0 0 1 
cosp 0 -—sing 
| 0 1 0 
sing 0 cosy 
cos sing 0 
E Ө соѕ Өө ^ 


0 0 1 


——| 


R,= 
For later purposes, let us write 


RSRQRQ,QRQR- P1 ә Рз (3.17) 
F1 F32 F33 
Up to now we have discussed rotations about vectors with origin O as their initial 
point. We now wish to drop the latter restriction and to perform a rotation about a 
vector v whose initial point can be any given point A(a;, a>, аз). To achieve this, we 
use vector v to compute matrix К of Eq. (3.17) in the same way as before. The 
following three steps are then taken: 


(1) Referring to Eq. (3.13), we perform a translation from the given point A to the 
origin O, using homogeneous co-ordinates and the following matrix: 


] 0 0 0 


0 1 0 0 

Te 
0 0 1 0 
=a =a —a, 1 


(2) We can now rotate about an axis passing through О as before, but we extend 
matrix R of Eq. (3.17) in a trivial way to be able to apply it to homogeneous 
co-ordinates: 


у 
T 
N 
mi 
دیا‎ 
= о о о 


EXERCISES 59 


(3) We apply a translation opposite to that in (1) using: 


1 0 0 0 
r-l9 1990 
“10 0 1 0 


The general rotation matrix is then 
Roen = TRT 
and it is used as follows: 


[x* y* gt 1] = [x y FA 1]Ксем 


EXERCISES 


3.1 Investigate the effect of choosing a = В = y = 1 in Section 3.6, Eq. (3.8) and 
write a program to perform the mapping from the quadrant to a given triangle 
for some points. Produce the triangle and the images of (1,1), (1,2), 
(1,3), . . . ; (2,1), (2, 2), (2,3),..., etc. in graphical form. 

3.2 Wirite a program for a general rotation. The input of the program will be: 


(1) The co-ordinates of two points A and B, determining the vector v — AB 
with initial point A. 

(2) The angle a. The rotation will then take place about vector v and through 
angle a. 

(3) A set of points to be rotated. 


The output of the program has to be the rotated points. Each point is 
expressed as a triple of rectangular co-ordinates. Since we have not yet dealt 
with graphical output of three-dimensional objects produce the output in 
numerical form. 

3.3 Find a polygon that will not be dealt with correctly by program POLY. TRIA, 
and develop a more general algorithm. 
Hint: If a triangle is cut off a polygon, none of the other vertices may lie inside 
that triangle. 


СНАРТЕК 4 


Perspective 


4.1 INTRODUCTION 


In Fig. 4.1 a two-dimensional representation of a cube is shown along with some 
auxiliary lines. In this picture lines such as AB and AD are not parallel to the lower 
and upper edges of the paper, so one could argue that they are not horizontal. 
However, they denote horizontal edges of the cube ABCDEFGH in three- 
dimensional space and we shall therefore briefly call them horizontal lines. For the 
same reason we shall say that two lines such as AB and DC are parallel, implicitly 
thinking in three dimensions. In this terminology, parallel horizontal lines meet in a 
so-called vanishing point. All these vanishing points lie on the same line, which is 
called the horizon. Note that the horizon and vanishing points are picture concepts, 
not really existing in three-dimensional space. For many centuries these concepts 
have been used by artists to draw realistic pictures of three-dimensional objects and 
it is customary to call such pictures perspective. 

The invention of photography offered a new (and easier) way of producing 
perspective pictures. There is a strong analogy between a camera used in 
photography and the human eye. Our eye is a very sophisticated instrument of 
which a camera is an imitation. In the following discussion the word eye may be 
replaced with camera if we wish to emphasize that two-dimensional hard copy is 
wanted. 

It is obvious that the picture will depend on the position of the eye. An important 
aspect is the distance between the eye and the object, since the ‘perspective effect’ 
will be inversely proportional to this distance. If the eye is close to the object we 
have a strong perspective effect as in Fig. 4.2(a). Here we can very clearly see that 
in the picture the extensions of parallel line segments will meet. On the other hand, 
if the eye is very far from the object (compared with the size of the object), parallel 
lines seem to be parallel in the picture. This is shown in Fig. 4.2(b). 

Besides the classical and the photographical methods there is a way of producing 
perspective pictures which is based on analytical geometry. We are by now very 
familiar with representing points in 2-space and 3-space by their co-ordinates (X, Y) 
and (х, у, 2), respectively. (Again we write ‘n-space’ as short-hand for ‘n- 
dimensional space’.) If we wish to produce a perspective drawing we are given a 
great many points P(x, y, z) of the object and we want their images P’(X, Y) in the 
picture. Thus all we need is a mapping from the so-called world co-ordinates 
(x, y, z) of a point P to the screen co-ordinates (X, Y) of its central projection Р’. 
We imagine a screen between the object and the eye E. For every point P of the 
object the line PE intersects the screen in point P'. It is convenient to perform this 
mapping in two stages. The first is called viewing transformation ; point P is left at its 
place but we change from world co-ordinates to so-called eye co-ordinates. The 
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Fig. 4.1. Vanishing points at the horizon 


СЕ‏ لل 


(a) (b) 


Fig. 4.2. (a) Eye close to object; (b) eye far from object 


second stage is called perspective transformation. This is a proper transformation 
from P to P', combined with a transition from the three-dimensional eye 
co-ordinates to the two-dimensional screen co-ordinates: 


World co-ordinates (Xw, Yw, Zw) 


^w 
Viewing аот 
Eye co-ordinates (Xe, Ye, Ze) 


Perspective vi) 


Screen co-ordinates (X, Y) 


4.2 THE VIEWING TRANSFORMATION 


To perform the viewing transformation the viewpoint E, being the position of the 
eye, and an object must be given. Let us require that the world-co-ordinate system 
be right-handed. It is convenient if it has its origin O lying more or less central in the 
object, since we then view the object from E to O. We shall assume that this is the 
case; in practice this might require a co-ordinate transformation consisting of 
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decrementing the original world co-ordinates by the co-ordinates of the central 
object point. We shall include this very simple co-ordinate transformation in our 
program, without writing it down in mathematical notation. 

Let the viewpoint E be given in spherical co-ordinates p, Ө, q, with respect to the 
world-co-ordinate system. Thus its world co-ordinates are 


Xg = p sin 9 cos Ө 
Yg = p sin ф sin Ө (4.1) 
ZE = p COS ф 


as shown in Fig. 4.3. The direction of vector EO (— —OE) is said to be the viewing 
direction. From our eye in E we can only see points within some cone whose axis is 
EO and whose top is E. If rectangular co-ordinates xg, yg, Zg were given, we could 
compute the spherical co-ordinates, as in Section 3.7. 

Our final objective will be to compute the screen co-ordinates X, Y, where we 
have an X-axis and a Y-axis lying in a screen between E and O and perpendicular to 
the viewing direction EO. This is why the eye-co-ordinate system, which we shall 
deal with first, will have its x.-axis and y.-axis perpendicular to EO, leaving the 
z,-axis in the direction of EO. The origin of the eye-co-ordinate system is viewpoint 
E, as shown in Fig. 4.4. Viewing from E to О the positive x.-axis points to the right 
and the positive y.-axis points upwards. These directions will enable us later to 
establish screen axes in the same directions. The z,-axis is chosen such that points 
with large positive z.-co-ordinates are far away. All this is logical and convenient, 
but its consequence is that the eye-co-ordinate system is left-handed. Although this 
might seem odd, it is quite usual in computer graphics and it will cause no problems 
at all. (Note that our world-co-ordinate system will always be right-handed.) 


Fig. 4.3. Spherical co-ordinates of viewpoint E 
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Fig. 4.4. Eye-co-ordinate system 


The viewing transformation can be written 
[X Je Az 185, we 2 ЫЎ (4.2) 


where V is the 4x4 viewing matrix. To find V we imagine the viewing 
transformation to be composed of four elementary transformations, for which the 
matrices can easily be written down. Matrix V will be the product of these four 
matrices. Each of the four transformations is in fact a change of co-ordinates and 
has therefore a matrix which is the inverse of the matrix that a similar point 
transformation would have. 


(1) Moving the origin from O to E 


We perform a translation of the co-ordinate system such that viewpoint E becomes 
the new origin. The matrix for this change of co-ordinates is: 


1 0 0 0 


к (4.3) 
0 0 1 O i 
1 


—Xp TYE TZE 


The new co-ordinate system is shown in Fig. 4.5. 


(2) Rotating the co-ordinate system about the z-axis 


Referring to Fig. 4.5 we now rotate the co-ordinate system about the z-axis through 
the angle л — Ө in the negative sense. This has the effect that the y-axis obtains the 
direction of the horizontal component of OE and that the x-axis becomes 
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Fig. 4.5. New axes after translation 


perpendicular to OE. The matrix for this change of co-ordinates is the same as for a 
rotation of points through the same angle in the positive sense. The 3 x 3 matrix for 
this rotation is: 


cos r — 0) sin (r —0) 0 
R,=| —sin (1л – Ө) cos (r 0) 0 


0 0 1 
sin 0 со$Ө 0 
—|-cos0 sin@ 0 (4.4) 
0 0 1 


The new position of the axes is shown in Fig. 4.6. 


| 
| 
| 
| 
ua 


N 
Fig. 4.6. New axes after rotation about z-axis 
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(3) Rotating the co-ordinate system about the x-axis 


Since the z-axis is to have the direction EO we now rotate the co-ordinate system 
about the x-axis through the angle л — ф in the positive direction. This corresponds 
to a rotation of points through the angle —(z — ф) = ф — x. Referring to Eq. (3.10) 
we have the following matrix: 


1 0 0 
R,=|0 cos(g—x) sin(g-—z) 
0 -sin(p—z) cos(qg — x) 
1 0 0 
=|0 -cosgp -—sing (4.5) 


0 sing —cos go 


The new axes are shown in Fig. 4.7. 


Fig. 4.7. New axes after rotation about x-axis 


(4) Changing the direction of the x-axis 


In Fig. 4.7 the y-axis and the z-axis have the right positions, but the x-axis is to point 


in the opposite direction. Thus we need the matrix for x' = —x, which is: 
-1 0 0 
M,=| 0 1 0 (4.6) 
0 0 1 


After this final transformation we have obtained the eye-co-ordinate system, already 
shown in Fig. 4.4. 
We can now compute the viewing matrix V as the matrix product 


V-TRZRiM; (4.7) 


where the notation R* is used for the 4 x 4 matrix obtained by extending a 3 x 3 


66 PERSPECTIVE 


matrix R with a fourth row and a fourth column containing the numbers 0, 0, 0, 1, in 
that order. Matrix multiplication is not commutative (in general, AB # BA), but it is 
associative, so we can write Eq. (4.7) as 
У = Т(К,К,М,,)* 

In this way we can deal with 3 x 3 matrices as long as possible. Our rather elaborate 
multiplication task is further relieved by introducing the abbreviations: 

cos p =a cos 0 = с 

sng-b sin 0 =d 
Thus à? + b? =1 and c? + d? =1. 

Using Eqs (4.1), we rewrite Eq. (4.3) as 


(4.8) 


1 0 0 0 

T= 0 1 0 0 
0 0 1 0 

—pbc -—pbd -pa 1 


and Eqs (4.4) and (4.5) as 


d c 0 
R,=|-c d 0 
0 0 1 
1 0 0 
R,=|0 -a -b 
0 b —a 
Hence 
d  —ac -—bc 
R,R,=| —c —ad -bd 
0 b —a 
Post-multiplying this matrix by M,, of Eq. (4.6) we obtain 
—d -—ac —bc 
R,R,M,,=| c -ad -bd 
b —a 
We then find the desired viewing matrix V as the product of 
1 0 0 0 
T- 0 1 0 0 
0 0 1 0 
—pbc -—pbd -pa 1 
and 
—d -—ac -—bc 0 
c =аа —bd 0 
R,R,M,,)* = 
( ya) 0 b на 0 
0 0 0 1 
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Hence 
-d -ac -bc 0 
c —ad —bd 0 
"uw ol ae 
Ug Uy Уаз 1 
where 


Ua = pbcd — pbcd = 0 

04 = pabc? + pabd? — pab 
= p(ab(c? + d?) — ab} 
— p(ab — ab) 
=0 

Vas = pb?c? + pb?d? + pa? 
= p{b*(c? + d?) + а?) 


= م‎ (b+ a?) 
=p 
Thus we have found 
—sin@ —соѕ ф соѕ0 -—singcos0 0 
y = cos Ө c g~sin@ -—singsinO 0 (4.9) 
0 sin Q —cos Q 0 
0 0 p 1 


We have now derived an important result. If we are given the spherical 
co-ordinates p, 0, q of viewpoint E we can compute the eye co-ordinates of a point 
from its world co-ordinates, using only Eqs (4.2) and (4.9). 

The viewing transformation, which we have now dealt with, is to be followed by 
the perspective transformation in the next section. However, we could also use the 
eye co-ordinates x. and y., simply ignoring z.. In that case we have a so-called 
orthographic projection. Every point P of the object is then projected into a point P' 
by drawing a line from P, perpendicular to the plane through the x-axis and the 
y-axis. It can also be regarded as the perspective picture we obtain if the viewpoint 
is infinitely far away. An example of such a picture is the cube in Fig. 4.2(b). 
Parallel lines remain parallel in pictures obtained by orthographic projection. In 
practice such pictures are quite usual. 

On the other hand, bringing some perspectivic effect into the picture will make it 
much more realistic. Our viewing transformation will therefore be followed by the 
perspective transformation, which will involve surprisingly little computation. 


4.3 THE PERSPECTIVE TRANSFORMATION 


The reader might have the impression that we are only half-way, and that this 
section will offer at least as much mathematics as Section 4.2. However, most of the 
work has been done. In this section the world co-ordinates will not be used. We 
shall therefore denote the eye co-ordinates by (x, y, z) instead of (Xe, Ye, Ze). 
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P(x,0, 2) 


Screen 


Fig. 4.8. Screen and eye co-ordinates 


In Fig. 4.8 we have chosen a point О, whose eye co-ordinates are (0, 0, d) for 
some positive number d. The plane z — d is the screen that we shall use. Thus the 
screen is the plane through O and perpendicular to the z-axis. The screen-co- 
ordinate system is then determined by prescribing its origin to be O and its X- and 
Y-axes to have the same directions as the x- and y-axes, respectively. For every 
object point P, the image point P' is the intersection of line PE and the screen. To 
keep Fig. 4.8 simple, we consider a point P whose y-co-ordinate is zero. However, 
the following equation to compute its screen co-ordinate X is also valid for other 
y-co-ordinates. In Fig. 4.8 the triangles EPR and ЕР'О' are similar. Hence 


PO PR 
EQ ER 
So we have 
ME. 
d z 
(4.10) 
X=d-~ 
2 
In the same way we can derive 
Y=d.? (4.11) 
2 


At the beginning of Section 4.2 we chose the origin О of the world-co-ordinate 
system to be a central point of the object. Since the z-axis of the eye-co-ordinate 
system is line EO, which intersects the screen in О, the origin О of the 
screen-co-ordinate system will be central in the image. If we require this origin to be 
in the bottom-left corner of the screen, and the dimensions of the screen are 2c, 
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horizontally and 2c; vertically, we have to replace Eqs (4.10) and (4.11) with 


E 94-2 (4.12) 
2 
T. 
ү=й-24а (4.13) 
2 


We still have to specify the distance d between viewpoint E and the screen. 
Roughly speaking, we have 


picture size object size 


ар 
which follows from the similarity of the triangles EP'O' апа ЕРО in Fig. 4.9. Thus 
we have 


be NAE He (4.14) 
object size 

This equation should be applied to both the horizontal and the vertical directions. It 
should be interpreted only as a means of obtaining an indication about an 
appropriate value for d, rather than as an exact prescription, since the three- 
dimensional object may have a complicated shape and it may not be clear how its 
size is to be measured. We then use a rough estimation of the object size, such as 
the maximum of its length, width and height. The picture size in Eq. (4.14) should 
be assumed to be somewhat smaller than the screen. More sophisticated methods to 
Obtain a desired picture size will be discussed in Section 4.5. 

The remaining part of this section is devoted to obtaining a better insight into 
concepts such as vanishing points and horizon. In Fig. 4.10 we have a viewpoint E 
and a screen ABCD. Again the viewing direction is from E to Q. Lines AF, BG, 
EH are parallel, and we consider them to be horizontal. We imagine a plane 


— р — 


Fig. 4.9. Picture size 
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Fig. 4.10. Vanishing point H 


Fig. 4.11. Vanishing point F of vertical lines 
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through the parallel lines EH and BG. This plane and the screen have BH as a line 
of intersection. Thus every point P on BG has its central projection P’ on BH, E 
being the centre of projection. If we let P move increasingly farther away from B, 
towards infinity, its projection P’ will approach H. This means that H is the 
vanishing point of the line through B and G. In terms of projective geometry H is 
the projection of the infinite point of BG. As discussed in Section 3.6, parallel lines 
are said to meet in an infinite point. In Fig. 4.10 the parallel lines BG and AF have 
the same infinite point, so H is also the projection of the infinite point of AF. If we 
take a line with a different direction, but also horizontal, such as line BF, this line 
also has its vanishing point on line CD. It is found as the intersection point of line 
CD and a line through E parallel to the given horizontal line. Line CD is the 
horizon. Every point J of the horizon is the vanishing point of all lines that are 
parallel to EJ. 

Not only horizontal lines have vanishing points. In Fig. 4.11 we have the vertical 
lines CA and DB. Point F is their vanishing point. It is the point where the vertical 
line through E meets the screen. Line segments CA and DB have projections CA’ 
and DB’, which are not parallel. 

We do not always appreciate a strong perspective effect for vertical lines but we 
often prefer pictures where vertical lines are represented by almost vertical lines. 
This is so because we are accustomed to horizontal or nearly horizontal viewing 
directions. Some people even get dizzy if they look a long way vertically! Artists 
produce ‘pseudo-perspective’ pictures where vertical lines are represented by exactly 
vertical lines, even if the viewing direction is not horizontal. In the latter case, such 
a picture differs from what we actually see but, curiously enough, it seems to be all 
right. An example is Fig. 4.12(b). In Section 4.6 we shall revisit this phenomenon, 
and show that our programs are capable of producing such pseudo-perspective 
pictures (see Fig. 4.22(a)). In general, however, it is a good idea to choose a 


(a) (b) 
Fig. 4.12. (a) Perspective; (b) pseudo-perspective 
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A B C D Е p G H 1 J 
Fig. 4.13. Ten cubes parallel to the screen 


viewpoint not too close to the object, especially not if the angle q indicated in Fig. 
4.3 differs much from 90°. This is a practical way avoiding a strong perspective effect 
for vertical lines. 

Another somewhat controversial subject is the representation of lines that are 
parallel to the screen. They will be represented by parallel lines in the picture. For 
example, consider Fig. 4.13, where we have ten cubes that are horizontally in line. 
The viewing direction is horizontal, so Фф = 90°. This picture has a paradoxical 
aspect, which has to do with the size of the cubes. Cube A is more remote than cube 
E and yet these two cubes are the same size in the picture. This seems to be 
incorrect since remote objects ought to be represented as being smaller than near 
ones. Yet it is a central projection, obtained by rigorous geometric rules deserving 
our confidence. The paradox is solved by observing that Fig. 4.13 is unrealistic with 
respect to the way we view objects. The construction of our eye is such that we only 
see points within a certain cone, whose axis is the viewing direction EO. An 
important aspect of this cone is the angle a. indicated іп Fig. 4.14. 

Eyes and cameras allow only a-values that do not exceed a certain aq. 
Computationally, however, we have not limited the actual angle а, indicated in Fig. 
4.14 and roughly satisfying 


0.5 - object size 
p 


so again a too-small choice of the viewing distance p is the source of the trouble. If 
we choose p large enough, the angle a will be sufficiently small to avoid the 
problems we discussed. If we are not yet satisfied there is another solution, namely 
to replace a flat screen with a part of a sphere whose centre is E. In this way we can 
allow large angles œ, but then we no longer have flat pictures. These spherical 
pictures themselves can be projected into a plane. Readers who are interested in 


tan œ = 


Fig. 4.14. Viewing cone 
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such unorthodox projection methods are referred to the work of the great graphical 
artist M. C. Escher (1972). We shall restrict ourselves to flat pictures and to rather 
small values of a. 


4.4 A PROGRAM TO DRAW A CUBE 


For many programming tasks we have to decide how general the program must be. 
There are two extreme possibilities. On the one hand, we can write a very special 
program, which does not read any input data and can only produce one result, say, 
some picture. On the other, we can develop a very general program, such as 
GENPLOT in Section 2.6, which can draw any picture provided that a file with 
input data be given. Between these extreme cases there are many possibilities for 
programs producing pictures of a certain type, after having read a limited amount of 
input data, sometimes called parameters. 

In this section we shall develop a program to draw a perspective picture of a cube 
whose edges have length 100. This completely arbitrary size is no restriction 
whatsoever with respect to the picture, since the following parameters can freely be 
chosen. 

(1) The three spherical co-ordinates (see Fig. 4.3): 
p, the viewing distance EO; 
0, an angle, measured horizontally from the x-axis; 
q, an angle, measured vertically from the z-axis. 

(2) The distance d between the screen and the viewpoint. 


Fig. 4.15. Cube and world co-ordinates 
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The origin O of the world-co-ordinate system is chosen at the centre of the cube, as 
in Fig. 4.15. The length of each edge is denoted by 2h, so h = 50. Then the vertices 
of the cube have the following co-ordinates: 


А( h, —h, —h) 
B( h, h,-h) 
C(-h, h,-h) 
D(-A, =h; =h) 
E( h,—h, h) 
F( h, h, h) 
G(-h, h, h) 
H(-h,-h, h) 


We shall draw a so-called wire-frame model, which means that we do not 
distinguish visible and hidden lines. In our program it is as if we are moving a pen in 
three-dimensional space. We shall use the functions mu(x, y, 2) and dw(x, y, z) ina 
similar manner as move(x, y) and draw(x, y) in two-dimensional space. In the 
functions mv and dw the function perspective will be used, which performs both the 
viewing and the perspective transformation. For the sake of efficiency all coefficients 
not depending on the particular points A to H are computed beforehand by the 
function coeff. Here is the complete program: 


/* CUBE: A wire frame model of a cube */ 

#include <math. h> 

float vii, vi2, v13, val, v22, v23, v32, v33, v43. 
screen dist, с1=4. 5, c2-3 5; 


main() 

{ float rho, theta, phi, h=50.0; 
printf("Viewing distance rho = ЕО: "); scanf("Zf", %тһо); 
printf("Give two angles, in degrees. Nn"); 
printf("Theta, measured horizontally from the x-axis: "); 
scanf("Xf", &thetadi 
printf("Phi, measured vertically from the z-axis: ")i 
scanf("Xf", &phidi 
printf("Distance from viewpoint to screen: "); 
scanf("Xf", £screen dist)i 
coeff(rho, theta, phi); 


initgr()i 

mv(h, -h, -h)i du(h, h; —h); /* draw AB */ 
du(-h, h, ~h); /* draw BC #/ 
du(-h, h: h); /* draw CG #/ 
dw(-h, -h, h); /* draw GH #/ 
du(h, -h, h); /* draw HE #/ 
du(h, -h, —һ); /* draw EA #/ 
mv(Xh, bh, —Һһ); du(h, h; h); /* draw BF #/ 
dw(-h, h, h); /# draw FG #/ 
mv(h, h; h); du(h, -h, h); /* draw FE */ 
mv(h, —h, -h); du(-h, —h, -h); /* draw AD #/ 
du(-h, he -h)i /* draw DC #/ 
mv(-h, -h, -h); du(-h, -h, Һһ); /* draw DH #/ 
endgr()i 
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coeff(rho, theta, phi) float rho, theta, phi; 

< float th, ph, costh, sinth, cosph, sinph, factor; 
factor=atan(1.0)/45. Oi 
/* Angles in radians: */ 


th-thetasxfactori 
costh=cos(th); 
cosph=cos(ph); 


ph=phitfactori 


sinth=sin(th)i 
sinph=sin(ph)i 


/* Elements of matrix V, see Eq. (4-9): */ 
vii--sinth; vi2--cosphscosth; vi3--sinph3tcosth; 
vei=costh; v22=-cosph#sinthi v23=-sinph#sinth; 
v32=sinphi v33--cosphi 
v43=rhoi 
> 


mv(x, у, 2) float х, у, 2; 

{ float X, Yi 
perspective(x, uy, z, ЯХ, &Y)i 
move(X, Y); 

} 


du(x, у, z) float x, у, Zi 

{ float X, Yi 
perspective(x, y, 2, &X, &Y)i 
drau(X, Y); 

} 


perspective(x, у, 2, pX, pY) float x, у, 2, #pX, #pY; 
€ float xe, ye, ze; 


/* Ege coordinates, computed as in Eq. (4-2): #/ 


xe = vil#¥x + valu; 

ye = vig*x + v22*y + v32*riij 

те = м13#х + v23"y + У93#:2 + v43; 

/* Screen coordinates, computed as */ 
/* in Eqs. (4-12) and (4-13): */ 
*pX = screen distsxe/ze + cli 

*pY = screen dist*ye/rie + c2; 


This program actually produced the cubes in Fig. 4.16. Note that in the second 
picture we view the cube from the bottom and that in the third picture the 
projections of sides and diagonals happen to lie on the same lines. These cubes can 
only be recognized with some difficulty. Things would be much better if we had 
drawn solid objects instead of wire-frame models. This will be the topic of Chapter 
5, where an algorithm is introduced to find out whether a line segment is visible 
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Fig. 4.16. Cubes seen from various viewpoints 
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from the viewpoint. However, if we take some care in choosing the viewpoint we 
can draw not too complicated solid objects by omitting hidden lines ourselves. For 
the cube in Fig. 4.15 this is done by choosing the viewpoint such that only the 
squares ABFE, BCGF, EFGH are visible and by omitting the edges AD, CD, DH. 
We then simply omit the three program lines immediately preceding endgr( ); in the 
program CUBE. The thus modified version of CUBE was actually used to produce 
Fig. 4.2(a), where p —200, 0 —30, ф —70, d —3, and Fig. 4.2(b), where р = 
200 000, 0 = 30, ф = 70, d = 3000. 


4.5 DRAWING WIRE-FRAME MODELS 


We shall now consider a program which can produce a perspective picture of any 
wire-frame model composed of finite straight line segments. By omitting certain line 
segments we can also produce pictures of simple solid objects, as we did for a cube 
in Section 4.4. As before, we shall use a viewpoint E and a ‘object point’ О, such 
that EO is the viewing direction. It would, however, not be convenient for the user 
to require this point O to be the origin of his world-co-ordinate system. We shall 
only require the user's z-axis to be vertical. The user will have to specify the 
co-ordinates Xo, yo, zo of object point O, expressed in his own co-ordinate system. 
When the user's co-ordinates are read, they are converted to 'internal world 
co-ordinates' x, y, z, with O as origin, according to: 


x = user's x-co-ordinate — Xo 

y — user's y-co-ordinate — yo 

z = user's z-co-ordinate — zo 
Then the user has to specify the spherical co-ordinates p, Ө, Ф, of viewpoint E, 
relative to the new origin O, as in Fig. 4.3. The distance d between viewpoint E and 
the screen must also be given. The position of the screen is then determined, since it 
is perpendicular to the viewing direction EO. To specify the line segments we again 


imagine movements in three-dimensional space that are associated with correspond- 
ing pen movements in the picture. For every move, a code 


0 — pen up 
1 = pen down 
is given, which was also used in the file A.SCRATCH in Section 2.6. To draw the 
two adjacent line segments PO and OR we have to specify three input lines of the 
following structure: 
Xp yp Zp 0 (move to P) 
Zo yo Zo 1 (draw PQ) 
XR Yr Zr 1 (draw OR) 
The right-hand side of each input line is comment, which actually may be present in 


the input file. 
Before presenting the program we shall consider an example of a complete set of 
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Fig. 4.17. View from positive x-axis 


input data. Let us consider Fig. 4.10, which is not completely trivial. If we draw it by 
hand it is a problem to determine the position of point M on line KL, such that in 
three-dimensional space EM is perpendicular to plane BCDA. Let us take K as the 
origin of the user's co-ordinate system. We let the positive x-axis of this system pass 
through B and the positive y-axis through L. Viewed from the positive x-axis the 
yKz-plane is shown in Fig. 4.17. Three-dimensionally, K is in the middle of segment 
AB, whose length is 10. We can now give the complete set of input data for the 
picture we wish to draw: 


0 2 2 (object point O) 
30 —15 75 15 (rho, theta, phi, d) 
5 0 0 0 (move to B) 
5 2 4 1 (draw BC) 
—5 2 4 1 (draw CD) 
—5 0 0 1 (draw DA) 
3 0 0 1 (draw AB) 
5 9 0 1 (draw ВС) 
=5 9 0 1 (draw GF) 
—5 0 0 1 (draw FA) 
0 -5 4 0 (move to E) 
0 2 4 1 (draw EH) 
0 0 0 1 (draw HK) 
0 9 0 1 (draw KL) 
0 —5 4 0 (move to E) 
0 3 0 1 (draw EM) 


78 PERSPECTIVE 


Fig. 4.18. Sample output of program GPERS 


We have found the value 15 for d by substituting p — 30, picture size — 7, object 
size — 14 in Eq. (4.14). Although the program still has to be presented, we show its 
output in Fig. 4.18, resulting from the above input data. 

The program will read from a file whose name is given as a program argument. 
This means that if the program name is GPERS and the file name is PLANES, we 
type 


GPERS PLANES 


instead of just GPERS in the command that starts the program execution. Programs 
with arguments usually have the line 


main(argc, argu) int argc; char *argv| |; 


at the beginning of the function main. The program name, here GPERS, is also 
regarded as an argument. The argument count argc is the number of arguments 
(argc — 2 in the example) and the argument vector argv contains pointers to the 
arguments; in the example we have 


argv[0] = "GPERS" 
argv[1] = "PLANES" 


The viewing and the perspective transformations are performed in the same way as 
in program CUBE in Section 4.4: 
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/* GPERS: General program for PERSpective */ 

*include <math. h> 

#include <stdio.h> 

float vii, vi2, vid, vai, v22, v23, v32, v3% v43, 
screen dist, с1=4. 5, c2=3. 5; 


main(argc, argv) int argc; char *argv[1i 

< float rho, theta, phi, х0, yO, 20, х, y, Zi 
int code 
FILE *f#fp; 
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if (argc != 2) X printf("No file name given"); exit (1)i Y 


if (#р=Ғореп(атдмС11, "т"), fp==NULL) 


€ printf ("File Xs does not exist", argvt11); exit(1); 


a; 


fscanf(fp, "Af Xf ZF", xO, &yO, &z0)i skipf(fp); 


fscanf(fp, "Xf Xf Xf Xf", rho, &theta, phi, £&screen dist); 


coeff(rho, theta, phi); 
initgr(); 
while (skipf(fp), fscanf(fp, "Xf Xf ZF hd", 


Yx: ец, &z, &code)) 


if (code) du(x-xO, y-yO, 2-20); else mv(x-xO, 
endgr()i 


skipf(fp) FILE #fpi < while (getc(fp) != ‘\n’);i } 


coeff(rho, theta, phi) float rho, theta, phi; 


€ float th, ph, costh, sinth, cosph, sinph, factori 


factor=atan(1.0)/45. Oi 

иж Angles in radians: #/ 

th-thetasfactor; ph=phi*factori 

costh-cos(th); sinth=sin(th)i 

cosph-cos(ph); sinph=sin(ph)i 

/* Elements of matrix V, see Eq. (4-9): #/ 

vli--sinth; vi2=-cosph#costhi vi3=-sinph#costh; 

vei=costh; ve2=-cosph#sinth; v23=-sinph#sinth; 
v32=sinphi v33=-cosph; 

v43=rho 


mv(x, y, 2) float x, Ч, 2; 

X float X, Yi 
perspective(x, у, 2, &X, ЬҮ); 
move(X, Y); 

+ 


du(x, y, 1) float х, у, 2; 

€ float X, Yi 
perspective(x, y, 2, EX, &Y); 
draw(X, Ү); 


1—70); 
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Perspective(x, y, 2, рХ, рү) float x, у, 1, *pX, #pY; 
€ float xe, ye, rei 


/* Eye coordinates, computed as in Eq. (4-2): +/ 
xe vli*x + v2l*yi 


ye vi2*x + v22*"y + v32*rii 

те у1З#х + v23*y  v33*i + v43; 

/* Screen coordinates, computed as */ 
/* in Eqs. (4-12) and (4-13): #/ 
#pX = screen dist*xe/ze + cli 

*pY = screen dist*ye/ze + c2i 


The function fscanf delivers the value 0 when a read attempt fails, that is, when 
the end of the file has been reached. Since 0 means false in a logical context, the 
while loop terminates correctly. 

The user of GPERS has to specify a value for the screen distance d. For the 
following reasons this is unsatisfactory: 


(1) The user is interested in the picture size rather than in the screen distance. 

(2) The ‘object size’ in Eq. (4.14) is a rather vague notion, so d cannot easily be 
determined in an exact way. 

(3) We sometimes wish to specify a rectangular viewport, which might be only a 
part of our screen. Then the picture should fit into this viewport. It would be a 
nuisance if we had to convert viewport dimensions into a screen distance. 


There are several ways of improving our program in this respect, namely: 


(1) Clipping the three-dimensional object against the pyramid whose top is the 
viewpoint E and whose base is a certain window, given in world co-ordinates. 
For this method, which we shall not employ, the reader is referred to Newman 
and Sproull (1979). 

(2) Clipping the picture two-dimensionally against a given viewport. We shall not 
use this method either. 

(3) Adjusting the size and the position automatically, such that the picture fits 
exactly in a given viewport. We used this method for two-dimensional objects in 
Section 2.6, where the program GENPLOT acted as a post-processor, and we 
can use this program again. Moreover, we can consult the last thirteen lines of 
the program CURVGEN in Section 2.6, where plot data are written to the file 
A.SCRATCH. 


We sball use method (3), using two files and two programs, as shown in Fig. 4.19. 
Program GPERSF will not read the screen distance d but simply use the value 
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program 
GPERSF 


file A. SCRATCH 


program 
GENPLOT 


graphical 
output 


Fig. 4.19. Flowchart for GPERSF and GENPLOT 


d — 1. This value is in fact irrelevant since the picture size computed by GENPLOT 
does not depend on d but on the given viewport. Recall that GENPLOT asks the 
user for the boundaries of the viewport. At this stage the reader might wonder if we 
can also get rid of the obligation to specify the object point O, since we have 
described O as a point central in the object and it seems that such a point can easily 


be computed automatically. 
that the user has to specify. 

The following program G 
screen distance d at the en 
omitted. 


However, there are two reasons to let O remain a point 
We shall discuss these reasons in Section 4.6. 

PERSF accepts the same input file as GPERS, but the 
d of the second line is ignored and may therefore be 


/* GPERSF: General program for PERSpective, producing the */ 
/# output File A. SCRATCH, to be read by GENPLOT */ 
include <math. h> 
#include <stdio.h> 
float vii, vi2, v13, v21, v22, v23, v32, v33, v43i 
main(argc, argv) int argc; char *argvlL]i 
X float rho, theta, phi, xO, yO, 20, x, y, 2; 
FILE sfpin, #fpout; 
struct { float X, Y; int code; } si 
if (argc != 2) X printf("No input file given"); exit(1); } 
if (fpin-fopen(argvt11," "т"), fpin==NULL) 
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€ printf("File Xs does not exist", argvt11)0; exit(1); 
+ 
fscanf(fpin, "Xf Xf Xf", 2х0, &yO, #20); skipf(fpin); 
fscanf(fpin, "Xf Xf Xf", rho, &theta, £phi)i 
coeff(rho, theta, phi); 
fpout=fopen("a. scratch", "u"); 
while (skipf(fpin), fscanf(fpin, "Xf Xf Xf Xd", 
£x, &y, &z, &s. code)) 
< perspect(x-xO, y-yO., 2-20, #5. Х, &s.Y); 
furite(&s, sizeof s, 1, fpout); 
} 
fclose(fpout)i 
} 


skipf(fpin) FILE #fpin; < while (getc(fpin) != ‘\n’)i } 


coeff(rho, theta, phi) float rho, theta, phi; 

€ float th, ph, costh, sinth, cosph, sinph, factor; 
factor=atan(1.0)/45. Oi 
/* Angles in radians: */ 
th-thetas*factor; ph=phi*factor; 
costh=cos(th)i sinth-sin(th); 
cosph-cos(ph); sinph=sin(ph); 
/* Elements of matrix V, see Eq. (4-9). #/ 
vii=-sinthi vi2=-cosph*costhi vi3--sinph3*costh; 
v2l-costhi v22=-cosph#sinthi v23=-sinph*sinthi 

уЗг=51прһ; v33=-cosph; 
v43=rhoi 
} 


perspect(x, yr 2, рХ, pY) float x, ys 2, #pX, +рү; 
{ float xe, ye, те; 


/* Ege coordinates, computed as in Eq. (4-2): */ 
xe = vil#¥x + veal*#yi 

ye = у12#х + v22*y + у32+2; 

ze = vi3tx + V23#y + vOG#z + v43; 


/* Screen coordinates, to be adjusted by GENPLOT +#/ 
*pX = хе/те; 
*pY уе/ е; 


The function perspect in this program has been derived from the function 
perspective of our previous programs CUBE and GPERS, by setting screen.dist — 1 
and cl = c2 = 0. (If other values had been chosen, the contents of file A.SCRA TCH 
would have been different but the picture produced by GENPLOT would have been 
the same!) 


4.6 VIEWING DIRECTION, INFINITY, VERTICAL LINES 


We shall use the program GPERSF to discuss some interesting new aspects. 
Suppose that an object is very large in the x-direction, such as the beam in Fig. 4.20, 
whose dimensions are 200 х 2 х 2. The interesting aspect of this example is the 
object point O that we must specify to determine the viewing direction EO. Up to 
now we have chosen O in the centre of the object. What we actually need, however, 
is a point O whose image point lies in the centre of the picture. Often this does not 
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Fig. 4.20. Long beam 


make much difference, but sometimes it does. The situation is illustrated in Fig. 
4.21, where in the picture point O' lies in the middle of Q'U', but the original point 
O does not lie in the middle of QU. If the middle of QU had been chosen for O, the 
viewing direction would have differed significantly from the current viewing 
direction. 

Returning to the beam in Fig. 4.20 it is clear that point O should not be chosen in 
the middle of the beam but much nearer to the eye. This is possible because in our 
programs the object point is given by the user rather than computed as the object 
centre. In the case of a landscape, finite line segments in the picture may result from 
infinite lines in reality. For such an infinite object it does not make sense to speak of 
an object centre, although there is a viewing direction EO, so it does make sense to 


Fig. 4.21. Object point not in object centre 
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define an object point O. The beam in Fig. 4.20 is not exactly a landscape, but its 
length suggests infinity; if this is not clear, this length should be given a much 
greater value than 200. In this beam we chose O(—15, 1, 1). Except for the letters 
P,Q,..., this picture was actually produced by the programs QPERSF and 
GENPLOT. The intersection of the beam with the plane x = —15 is denoted by 
ABC, and the intersection with x = —100 by A'B'C'. So the object point lies in 
ABC, but A'B'C' is in the middle of the beam. The complete set of input data for 


Fig. 4.20 is: 


—15 
30 


(a) 


Fig. 4.22. (a) Object point somewhat above the cube; (b) object point too far above the cube 
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(object point O) 
(rho, theta, phi) 
(move to Q) 
(draw QR) 
(draw RS) 
(draw SP) 
(draw PQ) 
(draw QU) 
(draw UV) 
(draw VW) 
(draw WS) 
(move to R) 
(draw RV) 
(move to A) 
(draw AB) 
(draw BC) 
(move to А”) 
(draw A'B’) 
(draw B'C') 


(b) 
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Another aspect deserving special attention is the representation of vertical lines, 
briefly discussed in Section 4.3. In Fig. 4.11 we saw the projected vertical lines meet 
in vanishing point F. We also compared the cube representations in Figs 4.12(a) and 
(b). If the object point O is placed in the centre of the cube our program cannot 
produce Fig. 4.12(b). However, we can place object point O above the cube such 
that the viewing direction EO is horizontal. The screen will then be vertical, since it 
is perpendicular to the viewing direction. This means that the vertical cube edges 
will be parallel to the screen. Their projections in the picture will therefore be 
exactly vertical. The curious aspect is that although the viewing direction is 
horizontal we view the upper plane of the cube from above. The picture that we 
obtain looks quite natural, unless we exaggerate and choose the viewpoint too high. 
Figures 4.22(a) and (b) show examples of both cases. 

The lines in Fig. 4.22(a) were drawn by GPERSF and GENPLOT, using the 
following input data: 


0 0 2 (point O) 

5 30 90 (rho, theta, phi) 

1 1 =1 0 (move to P) 
—1 l =] 1 (draw PQ) 
—1 1 1 1 (draw QU) 

1 1 1 1 (draw UT) 

1 Ll =] 1 (draw TP) 

i —L s] 1 (draw PR) 

L —1 1 1 (draw RW) 

1 1 1 1 (draw WT) 

1 —1 1 0 (move to W) 
—1 -1 1 1 (draw WV) 
—1 1 1 1 (draw VU) 

2 0 0 0 (x-axis) 

1 0 0 1 

0 2 0 0 (y-axis) 

0 1 0 1 

0 0 2 0 (z-axis) 

0 0 1 1 


Figure 4.22(b) was obtained by choosing point О much higher, namely at (0, 0, 6) 
instead of (0, 0, 2). This is no longer a realistic picture of a cube. Figure 4.22(a), on 
the other hand, is quite acceptable. Some will even prefer it to the cubes in Figs 
4.2(a) and 4.12(a), where vertical edges had a vanishing point. If a cube is 
perspectively drawn by hand it is customary to leave vertical lines vertical, as in Fig. 
4.22(a). The other representations are more difficult to produce manually. With the 
computer they are produced just as easily and the user can choose what he likes. 


REMARK 


The programs in this chapter are meant to explain principles rather than to be used 
in practice. In Chapters 5 and 6 we shall discuss more practical programs. 
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Fig. 4.23. Stairs 


Fig. 4.24. Pyramid 
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Fig. 4.25. Reduced cube 


EXERCISES 


In the following exercises hidden edges are to be omitted in the way we discussed at 

the end of Section 4.4. 

4.1 Write a program to draw the stairs of Fig. 4.23. The number of stairs (л) 
should be variable (n = 3 in Fig. 4.23). 

4.2 Wirite a program to draw a pyramid consisting of cubes, as in Fig. 4.24. The 
number 7 specifying the height of the pyramid in terms of cubes is to be read 
by the program (n —3 in Fig. 4.24). 

4.3 Write a program to draw the object of Fig. 4.25. This is a cube from which 
several small cubes have been removed. The dimensions of the small cubes are 
l/n times those of the original one. The number n is variable (n =3 in Fig. 
4.25). 


СНАРТЕК 5 


Hidden-line elimination 


5.4 CONSIDERATIONS OF EFFICIENCY 


Simple programs are not always faster than complex ones. For some problems 
simple algorithms are too slow and more sophisticated methods are required, 
sometimes resulting in quite complicated programs. Sorting is an example of such a 
problem. In Chapter 4 programs producing perspective drawings were rather simple 
and fast. Obviously the complexity of such programs will increase considerably if 
hidden line segments are to be eliminated automatically. Apart from the complexity 
that is inherent in the problem, additional complexity is inevitable if the programs 
have to be as fast as possible. Besides efficiency we can have requirements of 
generality, robustness and user's convenience as potential sources of program 
complexity. These four quality aspects deserve much attention, so the reader might 
fear that from now on we shall only deal with very complicated programs. It is, 
however, more attractive to begin with a rather straightforward program than with a 
fast but incomprehensible one. We shall therefore begin with a program that is of 
n^-time complexity. This means that computing time will be roughly proportional to 
n^, where n is the number of line segments or the number of polygons that have to 
be drawn. The program will be fast enough to draw simple objects, but it cannot be 
used practically for objects having many hundreds of faces. We shall store the data 
describing the object in arrays of fixed sizes. This will also limit the object's 
complexity. Apart from these limitations, the program will be quite general: in 
principle, all finite objects that have only a finite number of flat faces can be drawn. 
As to robustness and user's convenience, we shall not require that input errors will 
always result in clear messages or that the object can be specified in the most 
convenient way. Later we shall improve the program, especially with respect to 
efficiency. Hidden-line elimination is a notoriously time-consuming task and a real 
challenge for programming mathematicians. 


5.2 INPUT DATA AND INTERNAL REPRESENTATION 


An edge of an object may be entirely or partly hidden by one or more of that 
object's faces. Every edge is a (finite) line segment. An object may consist of several 
parts, which need not be connected. In Fig. 5.1 we have a tetrahedron and a cube, 
whose vertices are assigned numbers 0, 1, 2, .. . instead of letters A, B, C,.... The 
assigned numbers 12, 13, 14 will enable us to draw the positive co-ordinate axes. 
We intend to write a program that reads all data concerning the viewpoint and the 
object. For any (reasonable) viewpoint the graphical output will be the picture of 
the object. In contrast to Chapter 4, the object will be solid (and opaque). As in 
Section 4.5, the user has to supply two input lines specifying the rectangular 
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Fig. 5.1. Tetrahedron and cube 


co-ordinates of the object point О and the spherical co-ordinates р, 0, Фф of the 
viewpoint. In the example of Fig. 5.1 these input lines are 


2:5 1 1 (object centre O) 
8 20 70 (rho, theta, phi) 


We did not intend to pay much attention to the user's convenience but we shall 
introduce the concept of coded vertices: we shall need each vertex several times and 
it would be extremely annoying if three co-ordinates had to be given each time. In 
this example we specify the vertex numbers 0 to 11, each followed by the x-, y- and 
z-co-ordinates of the vertex. A separate line with the character # in its first position 
denotes the end of this part of the input: 


0500 (х0=5,у0=0,20=0) 8 2 0 2 
1320 (х1=3,у1=2,21=0) 9 2 2 2 
2 3 0 0 (ек) 10 0 2 2 
3 3 0 2 ng3s3 
4200 2700 
520 13 0 4 0 
6 0 2 0 14 0 0 3 
7000 # 


As in Section 4.5, the object point O is the origin of the world-co-ordinate system 
that is used internally. (Note the distinction between the letter О and the digit 0.) 
Since point О has co-ordinates (2.5, 1, 1), the co-ordinates of the points 0 to 14 are 
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decremented by these values. Thus point 9, for example, is internally given the 
reduced co-ordinates (—0.5, 1, 1) instead of the specified co-ordinates (2, 2, 2). 

Each face of an object can be any finite region of a plane, whose boundaries are 
line segments. A polygon is an example of such a region, but holes are allowed. 
However, to keep the program simple, the user has to decompose each region into 
triangles himself. We shall use these triangles to find out if they hide line segments. 
For reasons to be mentioned later we specify the numbers of the three vertices of 
each triangle in counter-clockwise order when viewed from outside the object. 
Every cube face is split up into a pair of triangles. Again the character # (on a new 
line) follows this piece of input data: 


1 3 0 (triangular faces of tetrahedron) 
1 2 3 
2 Q 3 
0 2. 1 
4 5 8 (front face of cube) 
8 5. 9 
5 6 9 (right-hand face) 
9 6 10 
8 9 10 (top face) 
8 10 11 
6 7 10 (back face) 
10 7 11 
7 4 8 (left-hand face) 
7 8 11 
6 5 4 (bottom face) 
6 4 7 
2 


Now every edge of the object has to be specified. Although these edges are 
already known as triangle sides they must be given again. First, not all triangle sides 
are object edges, and second, we appreciate the possibility of drawing additional 
line segments that do not belong to the solid object. To illustrate this we shall also 
draw parts of the positive co-ordinate axes, as in Fig. 5.1. Curiously enough, in the 
example this does not increase the number of input lines, since object edges that lie 
on the axes need not be specified again. Here are 17 input lines: 


4 12 (x-axis) 5 6 
7 13 (y-axis) 5 9 
7 14 (z-axis) 8 9 
0 1 (tetrahedron) 9 10 
1 3 10 11 
3. 10 8 11 
1 2 4 8 
2 3 6 10 
4 5 (cube) 


The program will read all input data from a file whose name is passed through the 
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program arguments argc and argu, as in Section 4.5. The co-ordinates of each vertex 
are stored in the array VERTEX, whose elements are structures containing the 
three fields x, y, z. The given user's co-ordinates are first reduced to co-ordinates of 
a system whose origin is the object point O. These internal world co-ordinates are in 
turn converted to eye co-ordinates, by means of the viewing transformation 
discussed in Section 4.2. It is these eye co-ordinates that are stored in array 
VERTEX: 


i VERTEX{i] 
SSS 


| x yz 
к 4o4d 


As we know, E is the origin of the eye-co-ordinate system and the viewing 
direction is the positive z-axis. Therefore all values VERTEX{i] - z must be positive. 
We shall check this in the program and stop if the check fails. 

A list of triangles is stored in the array TRIANGLE: 


j TRIANGLEY|j| 
| A B C a b e Ah 
L 4+ d 


Jl 1 4d 


Along with the vertex numbers of each triangle we store the coefficients a, b, c, h of 
the equation 


ax +by+cz=h (5.1) 


of the plane in which the triangle lies. We could calculate them each time they are 
needed but that would mean a waste of time, since we shall need them quite often. 
(We write h rather than d, since we have used d for the distance between viewpoint 
E and the screen.) We choose a, b and c such that 
a+b?+c=1 and hz0 
Then we can write Eq. (5.1) as the inner product 
n:x- 

where 

n-[a b c] 


х=[х y z] 
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The so-called normal vector n is the vector of unit length that is perpendicular to the 
plane of the triangle. For any point X in that plane, vector x = EX has the property 
that the inner product л = п. x is the distance between E and the plane. 

Inner products were introduced in Section 3.2. Note that we are using the 
eye-co-ordinate system, with E as origin and EO as the positive z-axis. A plane 
parallel to the screen has the equation z =h, which is a special case of Eq. (5.1), 
with a= b —0, c = 1. In the general case we compute a, b, c and h from the triangle 
vertices A, B, C. In Section 3.3 we saw that the plane through these points has the 
equation 


x y x l 
XA YA Za 1 
=0 
XB JB Zg 1 
Xo Ус & 1 
We can write this equation as 

Ya ZA 1 Xx Ze 1 XA Уд 1 ХА YA 24 
Ув Ze l|x—-|Xs 28 Ц y+ |хв ув 1 2= | хв ув Zn 
Ус 2с 1 Хс 2с 1 хс yc 1 Xc yc 2с 


In the program, a, b, c and h are therefore computed as follows: 


а= yA* (zB – 2С) — yB * (zA —zC) + yC * (zA — zB); 
b= —(xA*(zB —zC) - xB*(zA—zC)*- xC*(zA – 2В)); 
€ - xA*(yB — yC) - xB*(yA — yC) + xC*(yA – yB); 
h=xA*(yB*zC —yC*zB)— 

xB *(yA*zC —yC*zA) + 

xC *(yA*zB — yB * 2А); 
if (h 2 0) 
( r=sqrt(a*a+b*b+cx*c); 

a — alr; b —-blri; с= сін; h = hir; 


) else 

{... /* We can ignore triangle ABC for the reasons */ 
/* mentioned below. */ 

) 


We save a considerable amount of time by computing a, b, c and h only once, 
instead of each time that a line segment is investigated on visibility. 

If the program really had to be as simple as possible we would store all triangles 
that are read. However, triangles that are backfaces need not be stored but can be 
ignored. Consider triangle 123 in Fig. 5.1. In the picture this triangle is hidden by 
triangle 130. Since the latter triangle hides a part of line segment 45 we can ignore 
the fact that the former triangle does the same. A backface is hidden by visible 
faces. Although the backface can hide certain points from the eye, those visible 
faces also hide these points. This is why backfaces can be ignored. A simple way of 
determining whether a face is a backface is based on the orientation of the vertices. 
If we view the face 123 in Fig. 5.1 three-dimensionally from the outside (that is, 
from the negative x-axis) the order 123 in the input is counter-clockwise, since we 
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have requested this orientation in the input sequence. In the picture of Fig. 5.1, 
however, the order 123 is clockwise. This means that in the picture this face is 
viewed through the body of the object instead of from the outside. Thus 123 is a 
backface. It is a good exercise to apply this decision method to other triangles in 
Fig. 5.1. Referring to the end of Section 3.4, we find that we require the 
determinant 


Xa Ya 1 
D = Xp Ys 1 
Хе Yo 1 


where X4, Ya, XB, Yg, Xc, Үс are the screen co-ordinates of the vertices A, B, C of 
the triangle. Triangle ABC is a backface if D < 0. We do not actually compute D, 
since we can write 


XA/Za Yalza 1 XA YA ZA 
"” = |хв/2в ув/2в 1| = XB Ув Zp J льо) 
Xc/zc Yclzc 1 Xc Ус Zc 
= hl(zaZgZc) (5.2) 


Since z4, Zp, Zc are positive, D has the same sign as Л, computed in the part of the 
program above. If Л = 0, plane ABC passes through the origin E; but E is also our 
viewpoint, so in this case triangle ABC will hide no line segment and will therefore 
not be stored in array TRIANGLE. If h <0, determinant D is negative and triangle 
ABC is a backface. So we store triangle ABC in the array only if Л > 0, ог, in other 
words, if D >0. 

With the aid of vector products we could also have concluded directly from h <0 
that triangle ABC is a backface without using screen co-ordinates. We leave this as 
an exercise for the mathematically oriented reader. 


5.3 A HIDDEN-LINE ALGORITHM 


We shall now develop an important algorithm whose task is to draw a given line 
segment РО as far as it is visible. (We shall often briefly write ‘PQ’ to denote ‘line 
segment PQ’, to be distinguished from ‘line PQ’, which is the infinite line through 
the points P and Q.) For each triangle ABC of the list stored in array TRIANGLE 
we have to investigate whether it hides PO (or a part of it). 

As before, E is the viewpoint. We shall say that triangle ABC hides point R if line 
segment ER intersects triangle ABC in а point, interior to both the line segment and 
the triangle. The sides are not interior points of a triangle, and the endpoints are not 
internal points of a line segment. Thus R is not hidden by triangle ABC if R itself is 
an interior point of the triangle; neither are points on a side of the triangle hidden 
by the triangle. If triangle ABC does not hide point R we shall say that R is visible 
with respect to ABC. If at most a finite number of points of line segment PO are 
visible with respect to ABC, we shall say that ABC hides PO. For example, in Fig. 
5.1 triangle 130 hides line segment 12, although this triangle does not hide point 1. If 
a triangle does not hide PO we may not conclude that all points of PO are visible 
with respect to the triangle, since the triangle may partly hide PO. A triangle partly 
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hides PO if it hides infinitely many points of PO and if at the same time infinitely 
many points of РО are visible with respect to the triangle. 

If some triangle hides PO we need not take the remaining triangles of the list into 
consideration, but we can immediately conclude that РО is hidden and must not be 
drawn. 

We shall discuss the visibility of a line segment РО with respect to a triangle ABC 
in a more algorithmic way. We are given the eye co-ordinates xp, yp, Zp, Xo, etc. of 
the five points P, О, A, B, C and the coefficients a, b, c, h of the equation 
ax + by + cz =h that represents the plane through A, B, C. All information about 
triangle ABC is stored in the array element TRIANGLETj]. A key factor in our 
considerations will be the infinite pyramid, whose top is viewpoint E and whose 
planes pass through the sides AB, BC, CA of the triangle. This pyramid will briefly 
be referred to as ‘the pyramid’. In the same way the term ‘triangle’ refers to triangle 
ABC. 

Everything in the pyramid behind the triangle is invisible; all points in front of the 
triangle or outside the pyramid are visible (with respect to the triangle). In Fig. 5.2 
line segment РО intersects the pyramid in two points I and J. Subsegment IJ is 
invisible, PI and JQ are visible. (For the sake of brevity, we abbreviate ‘visible with 
respect to ABC’ to ‘visible’.) 

The complexity of our task lies in the great number of cases that we have to deal 
with. In the situation of Fig. 5.2 we can compute the position of I as follows. A 
vector representation of line РО is 


EP + Ar 
wherer-[rn n n]-PQ 
so rı = Хо – Хр 
№ = Уо ye 
F3 = 20 — 2р 


Thus point (x, у, 2) belongs to line РО if 


x = Xp + Ar 
y=ypt+ An (5.2) 
Z = Zp + Л» 


For À-values between 0 and 1, the point lies between P and О. 


Fig. 5.2. The pyramid 


A HIDDEN-LINE ALGORITHM 95 


The equation of plane EAB is 
x y z 
Xa YA Za| =0 
XB Ув 2в 
since it is an equation of a plane and it passes through the points E(0, 0,0), A, B. 
We write this equation as 
Cix + Gy + Ciz =0 (5.3) 
mum С, = уд2в — YBZA 
C; = XpZa — ХА2в 
C; = XAyB — Хвуд 
Substituting the right-hand sides of Eqs (5.2) into Eq. (5.3), we find: 


_ (Cixp + Сур + Cizp) 
(Cin + Gn + Gr) 


Using this À in Eq. (5.2) we obtain the co-ordinates of point I. The situation of Fig. 
5.2 occurs only if À lies between 0 and 1. For other values, point I does not lie on 
line segment PQ but on one of its extensions. We have not said that a value of A 
between 0 and 1 implies that PQ intersects the pyramid, but only the reverse. Figure 
5.3 shows a situation, viewed from E, where the value of А lies between 0 and 1, 
although PO does not intersect the pyramid. 

Returning to the three-dimensional world of Fig. 5.2 we imagine a plane through 
line РО and viewpoint E. We then want to know whether this plane passes through 
an interior point of AB. The calculations to find this out are similar to those that we 
used to find A. A vector representation of line AB is 


ЕА + uAB 


À = 


(5.4) 


Fig. 5.3. PQ outside pyramid 
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Thus points (x, y, z) belongs to line AB if 
X —XA t (Xp — XA) 
y = Ya + И(ув— ya) (5.5) 
Z* Za + u(zp — 24) 


For values of и between 0 and 1, the point lies between A and B. 
The equation of plane EPQ is 


x у z 
Xp yp Zp| =0 
Хо Yo Zo 
We write this equation as 
Kix + Ky + K3z =0 (5.6) 
where 
K = уь20 — YoZr 


К, = XoZp 7 Хр20 
K3 = Xpyo — хоур 
Combining Eqs (5.5 and (5.6) we obtain 


Kix, + Кул + K3za 


a чь — — Sel, 
K (xg — xA) + К(ув — Ya) + Кз(2в — Za) Ө.) 


p= 


Line segment PQ and pyramid plane EAB have a common point if and only if 
0xAx1 and Osu <1 


Up to now we have considered the intersection of PO with only one pyramid 
plane, namely EAB. We also have to take the planes EBC and ECA into 
consideration. Line segment РО may have zero, one or two intersections with the 
pyramid. The points P and O may lie inside, outside or on the pyramid; they may lie 
before, behind or on plane ABC. Such a check for one line segment and one 
triangle will involve a considerable amount of work. There are usually a great many 
line segments, and each of them has to be checked against a great number of 
triangles. These checks should therefore be done efficiently. (In Section 5.4 we shall 
improve efficiency with respect to the number of triangles to be matched with a 
given line segment, so we shall then reduce the number of checks. Each check for a 
given line segment and a given triangle, however, will then be carried out as here.) 
It is attractive to deal as soon as possible with cases that will often occur, especially 
if they do not require much computing time. We shall now list a sequence of tests 
that are performed in the given order. However, if a test succeeds, the remaining 
tests should be ignored: 


Test 1 (Fig. 5.4) 


If both points P and O lie before or in plane ABC (not behind it), PO is visible. This 
happens if 
n-EP<h and n-EQShA 
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Fig. 5.4. Points P and Q not behind plane ABC 


where n = [a b с]. Recall that a, b, с, h occurred in Eq. (5.1) and are available in 
array element TRIANGLETj]. Test 1 includes the important case that PQ is one of 
the sides of the triangle. 


Test 2 (Fig. 5.5) 


If (the infinite) line PO lies outside the pyramid, РО is visible. For this test we can 
substitute the co-ordinates of A, B and C in the left-hand side of Eq. (5.6), which 
denotes plane EPQ. If the three signs of the computed values (for A, B and C) are 
the same (all positive or all negative), points A, B, C lie on the same side of plane 
EPO. Then line PO lies outside the pyramid and is visible. We relax these sign 
conditions in the sense that one of the signs may be zero. We regard each sign as 
one of the numbers —1, 0, 1, and test whether the sum of the signs for A, B, C is 
one of the numbers —3, —2, 2, 3. Note that in Fig. 5.5(b) we could also have chosen 
point P coinciding with A. This, too, is an important special case. 


Q 
T" di C 
А 
(а) 


> 
[vo] 


B 


(b) 


Fig. 5.5. (a) |sign sum| = 3; (b) |sign sum| = 2 
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Test 3 


We now compute the intersections of line PO with the planes EAB, EBC, EAC. 
We compute 5 as in Eq. (5.4) for the intersection of line PO with plane EAB; we 
also compute и as in Eq. (5.7) for the intersection of line AB with plane EPO. (In 
both cases the line might be parallel to the plane. We shall then assign a very large 
value to A or u.) Equation (5.3) denotes plane EAB. If we substitute the 
co-ordinates of P and C into the left-hand side of this equation and the results for P 
and C have different signs, the points P and C lie on different sides of AB. We then 
say that P lies beyond AB. If a point lies beyond one of the sides of the triangle it 
lies outside the triangle. We store this information in the logical variable Poutside. 
The variable Qoutside has an analogous meaning. We shall use these variables in 
Tests 4 and 6. If P and O lie beyond the same side of the triangle, or if one of them 
lies beyond a side and the other precisely on that side, we say that PO lies outside 
the triangle. Line segment РО is then visible. Both situations are shown in Figs 
5.6(a) and (b). An important special case arises if Fig. 5.6(b) is somewhat modified, 
such that point О coincides with point A. Note that, unlike Fig. 5.5(b), the infinite 
line PQ in Fig. 5.6(b) may intersect BC between B and C. So Test 3 will succeed in 
cases where Test 2 has failed. 


C C 


Fig. 5.6. (a) Points P and Q beyond AB; (b) P beyond AB, Q on AB 


Most of the work has to be done three times, namely for each of the planes EAB, 
EBC and ECA, so we have three pairs (4,, ш) (L= 1, 2, 3). We then determine the 
minimum and the maximum of those values A, that satisfy the condition 


0xAx1 and 0<u<1 


and we call these values MIN and MAX, respectively. They are a by-product of this 
test and are defined only if PQ intersects the triangle, that is, if Test 3 fails and if 
Test 4, to be discussed now, also fails. 


Test 4 (Fig. 5.7) 


If neither Р nor О is outside the pyramid (and the previous tests have failed), PQ 
lies behind the triangle and is therefore invisible. 

This test also copes with a particularly tricky situation. Suppose, for example, that 
PQ lies on the pyramid behind AB. Since PQ then lies neither outside nor inside the 
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Fig. 5.7. PQ behind triangle ABC 


pyramid, it seems doubtful whether we may say that the triangle hides PQ. 
However, AB may not be an edge but a diagonal of a face. Then that face does hide 
РО, so РО must not be drawn. On the other hand, should AB be an edge of the 
object, then PQ need not be drawn, since in the picture PO would lie on AB and it 
does not make sense to draw coinciding line segments. 


Test 5 (Fig. 5.8) 


If a point I, where PO intersects the pyramid, lies in front of the triangle, PQ is 
visible. This test is based on the fact that the object is solid, so PQ cannot pass 
through an interior point of the triangle. To find such points I we use the A-values 
MIN and MAX, computed in Test 3. We can then simply test whether the inner 
product of the vectors EI and п = [a b c]isless than А. 


Q 


Fig. 5.8. Point of intersection in front of triangle ABC 


Test 6 


(This test is performed only if the previous tests have failed, that is, if PO intersects 
the pyramid behind the triangle.) 
If the A-values MIN and MAX yield the points of intersection I and J, as in Fig. 
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Fig. 5.9. Only one point of intersection 


5.2, then IJ is visible, and: 


If P is outside the pyramid or before plane ABC, PI is visible. 
If О is outside the pyramid or before plane ABC, JQ is visible. 


If P and О do not both lie outside the pyramid, then points I and J will coincide. 
This shown in Fig. 5.9. 

In the previous tests if PO was visible (with respect to the jth triangle), j was 
incremented by one. Line segment PO was then checked against the remaining 
triangles, if there were any, or drawn if there were not. In this test, however, there 
may be two sub-segments PI and JO that have to be checked against the remaining 
triangles. We can now recursively invoke our algorithm twice. For each recursive 
call, we specify the endpoints of the line segment and the number j +1 of the 
triangle where the remaining checks are to start. 


A note on numerical accuracy 


In complex computations it is wise to use type double instead of float. Since double 
stands for ‘double precision floating point’ we can informally speak of ‘floating 
point, if more specifically we mean double. Even if we use type double, real 
numbers are represented only with a finite precision. If the two numbers 5.843216 
and 5.843217 have been calculated in a complicated way we will probably have 
consideration for the difference in the last digit, and regard both numbers as equal. 
If we wish the computer to act similarly, we must specify this in a precise way. 
Ouantities of integer type are represented exactly. For floating-point types, 
however, this may not be the case. It is customary to choose a small positive value є 


(for example, ғ = 10~°) and replace the test x = = a (for floating-point values x and 
a) with 

fabs(x — a) <= epsilon, Or 

fabs(x — a) <= epsilon * a, Or 


fabs(x — a) <= epsilon + epsilon * a. 


In the first test, epsilon is the absolute error that we admit. This test may cause 
problems for large values of a and x unless we use a somewhat larger epsilon in 
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these cases. Since the magnitudes of x and a are usually unknown, we prefer a fixed 
value of epsilon, and consider replacing the first test with the second, where epsilon 
denotes the maximum relative error. However, if the value of a is (almost) zero the 
right-hand side of that inequality is (almost) zero, which is not exactly what we 
want. We shall therefore use the third test, where an absolute and relative error are 
combined. Obviously the third test reduces to the first if a — 0. 


Summarizing, we perform tests on floating-point values x and a as shown in the 
following table, where we use eps1 = epsilon + epsilon * a: 


Instead of We write 

х== fabs(x = a) < = ерѕ1 
х!=а fabs(x — a) > eps1 
x «a x «a — epsl 
х<=а x«-—a-epsl 
xa x >a + eps1 
x>=a X> = a — eps1 
х= =0 fabs(x) < = epsilon 
x!=0 fabs(x) > epsilon 
x <0 х < —epsilon 
x<=0 x < = epsilon 

х>0 x > epsilon 

x>=0 x >= —epsilon 


Note that the second half of this table can be derived from the first. In program 
HIDLIN below such tests will appear. 


/* HIDLIN: A simple program for hidden-line elimination. #/ 

/* The output of this program is the file А. SCRATCH, */ 

/* which is to be read by GENPLOT */ 

include “stdio һ2 

#include Zmath. h> 

#define NVERTEX 200 

#define NTRIANGLE 200 

int ntr=0; 

double vii, vi2, vi3, val, v22, v2% v32, v33. v4u. 
eps-le-9, meps--1e-5, oneminus-l1-i e-5, oneplus=i+1. e-5; 

FILE #fpout; 

struct € float X, Yi int code; У 5; 


struct < double x, ц, zi У VERTEX LNVERTEXJ; 
struct € int A, В, Ci double a, b. c, hi > 
TRIANGLE CNTRIANGLE]; 


main(argc, argv) int argc; char #argvl]i 
€ int i, ys А, В, С, Р, Gi 
double xO, yO, z0, rho, theta, phi, х, y 2, 
xà, ЧА, zA, xB, yB, zB, xC, YC, 2C, 
XA, ҮА, XB, YB, XC, ҮС, a, b, с, he rmi 
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FILE #fpin; 
char chi 


if (атдс!=2 ii (fpin-fopen(argvti11J, "r"))-zNULL) 
error("Input file not correctly specified"); 


fscanf(fpin, "AIF X1lf Alf", &xO, &yO, &z0); skipf(fpin); 
fscanf(fpin, "AIF ZIF X1f", &rho, &theta, &phi); 
coeff(rho, theta, phi); 
while (skipf(fpin), ch=getc(fpin), ch!=’#’) 
< ungetc(ch, fpin); 
fscanf(fpin, "74 X1f AIF X1lf", &i, Lx, By, &zDi 
if (i<O ii i>=NVERTEX) error("illegal vertex number"); 
viewing(x-xQ, y-yO, 2-20, 
&VERTEXCil. x, &VERTEXCil. у, &VERTEXE il. z); 
if (VERTEXCil. z<=eps) 
€ printf("Obyect point O and vertex Xd lie on ", i)i 
error("different sides of viewpoint Е. "); 
x 
F 


while (skipf(fpin), ch=getc(fpin), ch!=’#’) 

< ungetc(ch, fpin); 
fscanf(fpin, "Xd Xd 70", &A, &B, &C); 
xA=VERTEXCA]. xi yAZVERTEXEA1. у; zA=VERTEXCA]. zi 
хВ=УЕКТЕХЕВЈ. xi yB=VERTEXCB]. у; zBEVERTEXEB]1. zi 
xC-VERTEXECJ. xi yC=VERTEXCC]. yi zC=VERTEXCC]. z; 
a = yA * (zB-zC) = GB + (zA-zC) + ЧС + (zA-zB); 


b = -(xA * (zB-zC) - xB * (zA-12C) + xC + (zA-zB)); 
c = XA + (yB-yC) — xB * (yA-yC) + xC * (yA-yB); 
h = xà + (yB#zC = yC#zB) = 
xB # (yA#zC - yC#ZA) + 
xC # (yA#zB — yB*zA); 
if ¢(h>0) 
t if (ntr == NTRIANGLE) еттот ("Тоо many triangles"); 
T = sqrt(asatbsbtcstc); 
a = a/ri b = b/r; c = с/т; h= h/ri 
TRIANGLECntr].A = А; 
TRIANGLECntr1.B = Bi 
TRIANGLE(ntr1.C = C: 
TRIANGLEEntrl a = а, 
TRIANGLECntr1.b = b; 
TRIANGLE[ntr1.c = ci 
TRIANGLEEntr**J.h = hi 
У 
/* If h=0, plane ABC passes through Е and hides nothing. */ 
/* If h<O, triangle ABC is a backface. #/ 
/* In both cases ntr is not incremented and the triangle */ 
/* is not stored. */ 
F 
fpout-fopen("a.scratch", "о"); 


if (fpout==NULL) error("file a scratch cannot be opened"); 


while (skipf(fpin), fscanf(fpin, "7а 70", £P, &Q)) 
linesegment(VERTEXCPJ. x, VERTEXEPJ. y, VERTEXEP 1. 2, 
VERTEXCQ]. x, VERTEXEGJ. ч, VERTEXEGJ. z: О), 
fclose(fpout); 
} 


skipf(fpin) FILE #fpini { while (getc(fpin) '= ‘\n‘)i } 
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error(str) char *str; € printf("ZsNn", str); exit(10; + 
coeff(rho, theta, phi) double rho, theta, phi; 
€ double th, ph, costh, sinth, cosph, sinph, factor; 
factor=atan(i. 0) /45 0; 
/* Angles in radians: */ 
th=theta#factori ph=phi#factori 
costh=cos(th)i sinth=sin(th), 
cosph=cos(ph); sinph=sin(ph); 
/* Elements of matrix V, see Eq. (4-9): #/ 


vii=-sinth; 


v2l=costhi v22=-cosph#sinth; v23=-sinph#sinth; 
v32=sinphi v33=-cosph; 
v43=rhoi 
* 
viewing(x, y, 2, pre, руе, рте) 
double х, у, 2, #рхе, *pye, #pzei 
€ /* Eye coordinates, computed as in Eq. (4-2): */ 
*pxe = vil*x + v21*y; 
*pye = vli2*x + veeky + v32#%#z; 
*pze = vli3*x + у2З у + v33d*z + v43i 
> 
linesegment(xP, ЧР, zP, xG, YQ, 29, JO) 
double xP, ЧР, zP, xQ, Ча, 28; int JO; 
€ /* Line segment PG is to be drawn, as far as it is 
/* not hidden by the triangles JO, JO+1, ..., ntr 
int j=jO, worktodo=i, А, B, С, і, Pbeyond, Qbeyond, 
outside, Poutside, Qoutside, eA, eB, eC, sum, imin; 
double a, b, c, h, hP, hG, rl, r2, тз, 
хА, ЧА, zA, xB, yB, zB, xC, yC, 26, 
ЧА, dB, dC, MIN, MAX, lab, mu, 
xmin, ymin, тіп, xmax, утах, zmax: 
C1, C2, C3, Kl, K2, K3, denomi, denome, 
Cpos, Ppos, Gpos, aux, epsli 
while (y<ntr) 
< a=TRIANGLEC jJ. а; b=TRIANGLEL jJ. b; c=TRIANGLEL J1. c; 
h-TRIANGLEL jJ. h; epsi=eps+eps#h; 
/* Test 1: */ 
hP=a#xP+b#yP+c#2P; hG=a#xQ+b#yQG+c#zQ; 
if (hP<=h+epsl && hGz-h*tepsi) {j++; continue; } 
/# PG not behind ABC #/ 
/* Test 2: */ 
Ki=yP#zQ-yG#zP; K2-zP*xQ-zQ*xP; KG=xP#yQ-xQ#yP; 


A=TRIANGLEL jJl. Ai 
хА=УЕКТЕХГАЈ. xi 
xB=VERTEXCB]. xi 
xC=VERTEXCEC]. xi 
dA=K1#xA+tK2# yA+K3# 7 Ai 
dB=K1#xB+K2#*yBtK3#zBi 
dC=K1#xC+K2#yC+K3%zC; 


yA=VERTE 
yB=VERTE 


/* If dA, dB, dC have the 
lie at the same side of 
eA= dA>eps ? 1 dA<meps 
eB- dB>eps ? 1 dB<meps 
еС= dCreps ? 1 dC<meps 
sum = еА+еВ+еС; 


if (abs(sum)>=2) € Ј++; 


vi2=-cosph#costh; 


B=TRIANGLEL J1. В; 


yC=VERTEXCCI. y; 


vi3=-sinph#costh; 


C=TRIANGLEL J]. C; 
zA-VERTEXEA]1. 2; 
zB-VERTEXEB1. 2; 
z C-VERTEXEC2. z; 


ХГАЈ. yi 
XCBI. yi 


same sign, the vertices A, 
plane EPG. 
[o À OF 
2 oa s Ө» 
yt f OF 


continue; } 


*/ 
*/ 


B, 


*/ 
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/* If this test succeeds, the (infinite) line PG 
lies outside pyramid EABC 
(at most one common point) 
If the test fails, there is a point 
of intersection. */ 


/* Test 3: #/ 
Poutside-Goutside-O; MIN-1.; MAX=0. i 
for (1=0; i«3; i++) 
€ Cl=yA#zB-yB#zA;i C2-zAU5xB-zB*xA; C3=xA*#yB-xB#yAi 
/* C1 x + C2 y + СЗ z = О is plane EAB #/ 
Cpos=Cil#xC+C2#yC+C3#zC; 
Ppos=C1#xP+C2#yP+C3#zP; 
Qpos=C1#xQ+C2#yQt+C3#72Qi 
denomi=Gpos-Ppos; 
if (Cpos>eps) 
X Pbeyond- Pposcmeps; Gbeyond- Gpos<meps; 
outside= Pbeyond && Qpos<=eps if 
Gbeyond && Ppos<=eps; 
* else if (Cpos<meps) 
< Pbeyond= Ppos>epsi Qbeyond= Gpos>eps; 
outside= Pbeyond && Gpos>=meps ii 
Qbeyond && Ppos>=meps; 
+} else outside=1; 
if (outside) breaki 


lab= fabs(denoml)<=eps ? 1.e7 : -Ppos/denomli 
/* lab indicates where PG meets plane EAB #/ 
Poutside i= Pbeyond; 

Goutside i= Qbeyondi 

denom2=dB-dA; 

mu= fabs(denom2)<=eps ? 1. е7 : -dA/denome; 


/* mu tells where AB meets plane EPG #/ 
if (mu>=meps && mu<=oneplus && 
lab>=meps && lab<=oneplus) 
€ if (lab<MIN) MIN-labi 
if (lab>MAX) MAX-lab; 
F 
aux=xAi xA-xB; xB=xC; xC=aux; 
aux=yAi yA=yBi yB=yCi yC=auxi 
aux=1Ai zA=zBi zB=zCi zC-auxi 
aux=dAi dA=dBi dB=dCi dC=aux; 
+ 
if (outside) {j++; continue; } 


/* Test 4: */ 


if (!(Poutside |! Goutside)) 
< worktodo=0; break; /# PG invisible */ 
$ 


/* Test 5: */ 
ri-xG-xP; r2-yG-yP; r3-zGQG-zP; 
xmin=xP+MIN#r1; ymin=yP+MIN#r2; zmin=zP+MIN#r3i 
if (atxmintb#ymintc#zmin«h-epsi) < j++; continue; У 
xmax=xP+MAX#r1; ymax=yP+MAX#r2; zmax=zP+MAX#r3i 
if (a*xmax+b#ymaxtc#zmax<h-epsl) < j++; continue; > 


/* If this test succeeds, an intersection of PG 
and the pyramid lies in front of plane ABC. */ 
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/* Test 6: */ 


if (Poutside ii hPzh-epsi) 
linesegment(xP, yP, zP, xmin, 
if (Goutside ii hGcXh-epsi) 


linesegment(xG, Ча, zQ, xmax, 
worktodo=0; break; 
Y 


if (worktodo) 

€ s. X=xP/zPi 5. Ү=уР/тР; s.code=0; 
furite(&s, sizeof s, 1, fpout); 
S. X-xG/z2G0; s. Y-yG/:zGi s.code-l; 
furite(&s, sizeof s, 1, fpout); 


ymin, 


ymax, 


zmin, 


zmax, 


gti 


gti): 
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The input data of Section 5.2 was read by this program. The results written in the 


file A.SCRATCH were read by program GENPLOT, which produced Fig. 5.10. 


Fig. 5.10. Output of program HIDLIN 


5.4 POLYGONS AND PIXELS 


We shall now improve program HIDLIN of Section 5.3 with respect to the user's 
convenience and to efficiency. We shall enable the user to specify object faces as 
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Fig. 5.11. Solid letter A 


polygons' and the program itself will decompose these polygons into triangles. Line 
segments are then no longer to be given separately if they are edges of polygons, 
which is usually the case. 

Let us, for example, consider the solid letter A of Fig. 5.11. The front face and 
the back face are not proper polygons since they have a triangular hole. АП other 
faces are rectangles. The vertices in the front face have numbers 0 to 9. Increasing 
each number by 10 yields the number of the corresponding vertex in the backface. 
In the input file for our new program we can specify the front face as follows. 


0123456-9879 —6 


The minus signs in the sub-sequences 6 —9 and 9 —6 mean that these sub- 
sequences do not denote line segments to be drawn. They serve merely as a means 
of defining the region of the front face. The dashed line 6 9 in Fig. 5.11 can be 
regarded as two coinciding edges of a polygon. We obtain the front face in Fig. 5.11 
if the distance d in the polygon of Fig. 5.12 becomes increasingly smaller. Instead of 
6 9, other pairs of points could have been chosen to connect the outer and the inner 
boundaries. (However, we shall not write —0 since this does not represent а 
negative number.) The vertices must be listed in counter-clockwise order, with the 
exception of the inner boundary, where they are listed clockwise. The general rule is 


' See the footnote in Section 3-5. 
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Fig. 5.12. A polygon 


as follows. If we track the edges from vertex to vertex in the listed order, each time 
facing the next vertex, the region we are defining is to be on our left-hand side. This 
is why —9 in the above input sequence is followed by 8. If we move from vertex 9 to 
vertex 8 in Fig. 5.11, facing points 8, the ‘material’ of the letter A is on our left-hand 
side. Note that the input sequence is interpreted cyclically: the final vertex (6) is 
connected to the first vertex (0). 

We shall require the first three vertices to enclose a convex angle. For example, 
we see in Fig. 5.11 that 


123 
at the beginning of the sequence would have been illegal, but 
601 


would have been acceptable. 

The final vertex of every sequence must immediately be followed by the character 
#. This convention enables us to use several input lines for one sequence, so we can 
have polygons with a great many vertices. If a sequence contains only two vertices it 
is interpreted as a loose line segment. We use this facility, for example, to draw 
co-ordinate axes, as in Fig. 5.10. 

As before, the file will begin with the co-ordinates 


Xyz 
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of the central object point O, which will be used as a new origin. The spherical 
co-ordinates p, Ө, ф of viewpoint E, relative to O (see Fig. 4.3), will also be needed 
again, but instead of reading them from a file the program will now ask the user for 
these three numbers. In this way we can execute the program several times with 
different viewpoints without changing the file. Line EO is the viewing direction. The 
number of each vertex and its x-, y- and z-co-ordinates are again read from the file. 
Finally the vertex numbers of the faces are listed as we have discussed. Remember 
that we must view a face from outside the object to determine what is counter- 
clockwise. We shall adopt the convention that comment may appear anywhere 
outside numbers and that all comment be in parentheses. These parentheses must 
not be nested. To be able to determine where the last part of the input data begins it 
shall be preceded by the word Faces. (In fact we shall test only the first letter, 
capital F, and skip the remaining letters, so the word Facets could also be used.) 
The complete set of input data to draw the letter A of Fig. 5.11 is listed below. 


0 0 30 (co-ordinates of point O; 
the viewing direction is EO) 


0 0 —30 0 (vertex 0) 
10 —10 —30 0 (vertex 10) 
1 0 —20 0 (vertex 1) 
11 —10 -20 0 (vertex 11) 
2 0 —16 8 (vertex 2) 
12 —10 -16 8 (vertex 12) 
3 0 16 8 (vertex 3) 
13 —10 16 8 (уегїех 13) 
4 0 20 0 (vertex 4) 
14 —10 20 0 (vertex 14) 
5 0 30 0 (vertex 5) 
15 —10 30 0 (vertex 15) 
6 0 0 60 (vertex 6) 
16 —10 0 60 (vertex 16) 


7 0 —12 16 (vertex 7) 
17 —10 -12 16 (vertex 17) 
8 0 12 16 (vertex 8) 
18 —10 12 16 (vertex 18) 
9 0 0 40 (vertex 9) 
19 —10 0 40 (vertex 19) 


Faces: 


0 1 2 3 4 5 6-9 8 7 9 —6# 
10 16 —19 17 18 19 —16 15 14 13 12 11# 
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The program that we shall develop will be much more convenient than program 
HIDLIN in Section 5.3 since instead of triangles and line segments we can now 
directly specify object faces with complicated shapes. However, this is not the only 
improvement. In program HIDLIN we had a set of line segments and a set of 
triangles. Each line segment was matched with each triangle. Thus if the object 
increases in complexity in the sense that both sets become twice as large, computing 
time will increase by a factor 4. We say that HIDLIN has a quadratic time 
complexity. 

For complex objects our new program will be much faster than HIDLIN. As 
before, we shall have a set of line segments and a set of triangles, but a given line 
segment will now be matched with only a sub-set of the set of triangles. For 
example, a line segment PO that lies at the top left-hand part of the screen will not 
be matched with a triangle ABC in the bottom right-hand corner. To achieve this, 
we divide the screen into Nscreen X Nscreen rectangles of equal size, where Nscreen 
is some positive integer. In Fig. 5.13 we have chosen Nscreen = 8. 

Such an elementary rectangle is called a pixel, which stands for picture element. 
Pixels are usually associated with raster-scan displays. To have an acceptable 
resolution, with such hardware many hundreds of pixels in both the horizontal and 
the vertical directions are required, since each pixel is then entirely filled with one 
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Fig. 5.13. Device-independent pixels 
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colour. Up to now, our approach has been device-independent. Perhaps contrary to 
what the reader will expect after the introduction of pixels, we shall stick to this 
principle! We are in fact using device-independent picture elements, so we could 
have invented a new name, such as dipels. Instead we shall conform to the usual 
terminology and use the term pixel. Note, however, that our pixels have nothing to 
do with resolution or with hardware. The choice of Nscreen (8 in Fig. 5.13) may 
influence computation time, but it will not influence the resulting picture. Good 
performance figures have been obtained by the value Nscreen — 30, which means 
that there are 30 x 30 — 900 pixels. 

Our improvement of efficiency is based on the idea of setting up a list of triangles 
for each pixel. Such a pixel list will contain only those triangles that (partly or 
entirely) cover the pixel. Thus in Fig. 5.13 only the lists of the shaded pixels in 
screen columns 5, 6, 7 will contain triangle ABC. As in program HIDLIN, there 
will be one general array TRIANGLE, whose elements are structures containing the 
vertex numbers A, B, C and the coefficients a, b, c, А of the equation of plane 
ABC. Recall that the eye co-ordinates of vertex A are stored in array element 
VERTEX[A], etc. Depending on the context, ‘triangle ABC’ means either the 
original triangle in three-dimensional space or its projection onto the screen. We 
shall say, for example, that triangle side AB in Fig. 5.13 lies in pixel columns 5, 6, 7; 
this will somewhat simplify our discussion, since otherwise we would have to talk 
about the central projection A'B' of the original side AB. In the pixel lists we store 
only subscript values j, representing the triangle whose relevant characteristics can 
be found in TRIANGLE{j]. We shall say that the triangles that partly or entirely 
cover a pixel are associated with that pixel. When all pixel lists have been 
established we can start drawing line segments. For each line segment PO we begin 
with building the set of all triangles associated with pixels that contain points of PO. 
In Fig. 5.13 such pixels are shaded in columns 0, 1, 2. We shall then match PO only 
with the triangles in this set. At first sight it seems that the realization of this rather 
simple idea will require a considerable amount of additional work, since first, all 
triangles have to be stored in pixel lists and second, the visibility checks for each line 
segment have to be preceded by constructing the set of triangles that might hide it. 
Note, however, that the amount of work for these two actions depends only linearly 
on the number of triangles and line segments. Thus we can expect linear time 
complexity for this method. 

In Section 5.3 the viewport co-ordinates were not computed by HIDLIN but by 
the post-processor GENPLOT. The latter program first determines the minimum 
and the maximum values of the given co-ordinates (read from the file 
A.SCRATCH ) and uses these to compute the scale factor f and the offsets c1 and 
c2, as discussed in Section 2.6. We now prefer to incorporate this work into our new 
program HIDLINPIX. We need the minimum and maximum co-ordinate values 
anyhow, because otherwise we would not be able to associate screen points with 
pixels. It would be a waste of time to do this work twice. Internally we shall use 
screen co-ordinates X and Y, computed according to the ‘clean’ perspective 
transformations of Eqs (4.10) and (4.11), which means that this screen is at a 
distance 1 from viewpoint E. We compute these screen co-ordinates as soon as 
possible to determine their minimum and maximum values Xmin, Xmax, Ymin, 
Ymax. This screen should not be confused with the viewport where the picture will 


POLYGONS AND PIXELS 111 


eventually appear. The correspondence between screen and viewport follows from 
the viewport boundaries Xvp.min, Xup-max, Yup-min, Yvp.max, that the user is 
asked for when the program is executed. The following program fragment is closely 
related to Section 2.6 and shows what is actually computed: 


Xrange — Xmax — Xmin; 

Yrange = Ymax — Ymin, 

Xvp.range = Xvp.max — Xvup-min; 
Yup-range = Yup-max — Yup-min, 

fx = Xup-range/Xrange; 

fy = Yup-range/Yrange; 

f — (x « fy ? fx : fy); 

Xcentre = 0.5 * (Xmin + Xmax); 

Ycentre = 0.5 * (Ymin + Ymax); 
Xvp.centre = 0.5* (Xup-min + Xvp.max); 
Yup.centre = 0.5 * (Yup-min + Yvp.max); 
c1 = Xup-centre — f * X-centre; 

c2 = Yup-centre — f * Y centre; 


The viewport co-ordinates Xvp and Yup can then be computed from the screen 
co-ordinates X and Y as follows: 


Хор =f*X +с1; 
Yup =f *Y +c2; 


It is interesting to compare this with Section 4.3, which reveals that scale factor f 
is nothing else but the distance d between the viewpoint E and the plane where we 
can imagine the viewport. Thus the screen and the viewport are two parallel planes 
(both perpendicular to the viewing direction EO), the former at a distance 1 and the 
latter at a distance d = f from viewpoint E. 

We shall use the viewport co-ordinates only in the actual ‘plotting commands’, 
that is, in the calls of the functions move and draw. Internally we use screen 
co-ordinates X and Y, lying between the computed boundaries Xmin, Xmax, Ymin, 
Ymax. It is this screen that is divided into Nscreen X Nscreen pixels. We shall now 
discuss the relationship between x and ipix, the column number of the pixel. The 
relationship between Y and the row number jpix will be analogous. As Fig. 5.13 
shows, we have 


Xmin < X < Xmax 
0 € ipix < Nscreen — 1 
Theoretically, the horizontal dimension of a pixel is 
deltaX — (Xmax — Xmin)/Nscreen 


(We shall presently see that a slight correction to this equation is practical.) Then 
ipix is obtained by truncating the value of 


(X — Xmin)/deltaX 


to an integer. This computation has the unpleasant consequence that the value 
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X — Xmax will yield ipix — Nscreen, instead of the largest value Nscreen — 1 that is 
permitted. One solution would be to extend the screen with a column (and a row) 
which then would only be used for points on the line Х = Xmax (or Ү = Ymax). For 
reasons of economy we shall not use this solution, but insist on having only Nscreen 
pixel columns, ranging from 0 to Nscreen — 1. This can easily be accomplished by 
slightly increasing deltaX. We shall multiply it by 1+ epsilon, where, as usual, 
epsilon has a small positive value, say 107°. Thus, as soon as Xmin, Xmax, Ymin, 
Ymax are determined, we compute: 


deltaX = (1 + epsilon) * Xrange/Nscreen ; 
deltaY = (1 + epsilon) * Yrange/Nscreen; 


Then any time we want the pixel associated with point (X, Y) its position is obtained 
by: 


ipix = (X — Xmin)/deltax; 
jpix = (Y — Ymin)/deltaY ; 


Note that truncation is performed implicitly in the C language whenever a 
floating-point value is assigned to an integer variable. Thus, declaring ipix and jpix 
of type int is all we have to do. 

We shall also require the reverse operation: 


X = Xmin + ipix * deltaX; 


We thus find the X that corresponds to the left-hand boundary of pixel column ipix. 
The screen point (X, Y) in the centre of pixel (ipix, jpix) is computed as follows: 


X = Xmin + (ipix + 0.5) * deltaX; 
Y = Ymin + (jpix + 0.5) * deltaY; 


We shall now see how pixels can be used to store information about the triangles 
they are associated with. We use linear lists of nodes, and each node contains a 
triangle number jtr and a pointer to the next node, if there is one. If not, the pointer 
field will have the value NULL. The start pointer of each list is located in the 
two-dimensional array SCREEN. Pixel (ipix,jpix) corresponds to 
SCREEN|ipix J[jpix ]. 

For the sake of efficiency, some of the triangles associated with a pixel will be 
given special treatment. These special triangles are those that completely cover the 
pixel. (A triangle covers a pixel completely if the entire pixel lies within the 
triangle.) For a given pixel (ipix, jpix) all these special triangles can be ignored 
except the one that is nearest to the viewpoint! Here the terms ‘near’ and ‘distance 
of the triangle’ refer to the point in three-dimensional space where the line through 
E and the centre of the pixel meets the triangle. If there are triangles that 
completely cover the pixel, the number identifying the nearest of them will also be 
stored in SCREEN([ipix |[уріх], along with the distance of this triangle. The fields for 
these two new items are called tr.cov and tr.dist. Figure 5.14 shows a typical set of 
triangles associated with a pixel. 

During the process of storing triangles, the pixel can have the data structure 
shown in Fig. 5.15. Only triangles 18 and 23 cover the pixel completely. Since 
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18, near by 
(dist. = 25.7) 


23, far away 
(dist. = 80.3) 


Fig. 5.14. Pixel and associated triangles 


triangle 18 is the one that is nearer to the eye, its number and its distance are stored 
in the SCREEN element. Triangles 3 and 5 do not completely cover the pixel, 
therefore their distance does not matter at this stage. During the process of storing 
triangles, the nearest covering triangle, such as 18 in the example, is not 
immediately stored in the linear list since there might follow another triangle that 
covers the pixel and is nearer. At the end of this process the nearest covering 
triangle is added to the list. So finally the list of each pixel will contain all triangles 
that are associated with the pixel. 

We now have to cope with an interesting programming problem. For a given 
triangle ABC, with screen co-ordinates XA, YA, XB, YB, XC, YC, we have to find 
out which pixels are associated with it and which of them, if any, are completely 
covered by it. This problem is far from trivial, because there are so many cases to 
consider. The first sub-problem is to determine for each of the three triangle sides 
whether it is a boundary at the top or at the bottom of the triangle. We assign the 
numbers 0, 1, 2 to the sides AB, AC, BC, respectively, as in Fig. 5.16. Here side 0 
happens to be at the top and sides 1 and 2 at the bottom. Our intention is to code 
this as follows in array topcode, which has three elements: 


topcode[0] = 1 
topcode|1] = 0 
topcode|2] = 0 


So 1 and 0 are the codes for top and bottom, respectively. To determine these values 


SCREEN[ipi x][jpix] linear list 
trcov tr.dist start Jtr next jtr next 
s T H——B = 


Fig. 5.15. Data structure for pixel 
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А 


B 2 C 
Fig. 5.16. Side numbering 


analytically we use a point V that lies far below the screen. Let us denote the screen 
co-ordinates of V by (0, M), where, for example, M = —10°. Then topcode[0] = 1 if 
and only if both points C and V lie on the same side of line AB. The latter can easily 
be checked by substituting the co-ordinates of C and V into the left-hand side of the 


following equation of line AB: 


ХА Ya 1 
Хв Ув 1| =0 
zo у 1 


If C and V yield the same sign in this substitution they lie on the same side of AB, 
which implies that AB is a top side. Thus in the C language we can briefly write 


topcode[0] = (D * DAB > 0); 
where D and DAB denote the determinants 


Xa YA 
D-|xs ув 


Рв = |Хв ув 


= 
[e] 
ме 
Oo 
mM mom نم نم‎ pa 


Analogous to this, we compute the determinants 
XA Ya 1 

Dac-|0 M 1 

Xc yc 1 

0 M 1 

ЮОвс= |хв ув 1 

1 


Xe Ye 
and compute 


topcode|1] = (D * DAC >0); 
topcode[2] = (D * РВС > 0); 
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Our next step is to determine the left-hand endpoint (X/eft[/], Yleft[/]) and the 
right-hand endpoint (Xright[/], Yright[{]) of each side / (I = 0, 1, 2). For example, in 
Fig. 5.16 we have Xleft[0] = XB, Yleft[0] = YB, etc. 

We also determine the numbers ipixmin and ipixmax of the pixel columns where 
the extreme left and right vertices of triangle ABC lie. So for triangle ABC we have 
to update the pixel lists only for pixels in the columns ipixmin, . . . , ipixmax. For 
each column ipix in this range we introduce the row boundaries LOWER[ipix] and 
UPPER{ipix]. All pixels LOWER([ipix], ..., UPPER([ipix] are associated with the 
triangle. A (possibly empty) sub-range of this is LOW[ipix]+1,..., UP[ipix] — 1. 
All pixels in this sub-range are completely covered by the triangle. Figure 5.17 
shows the /th side of triangle ABC. Let us assume that this side is a bottom 
boundary, so topcode[/] = 0. Then this triangle side will contribute to the values of 
LOWER[ipi| and ГОМ [іріх] for the values of іріх in the range 
ipixleft,..., ipixright, where the integers ipixleft and ipixright are truncated 
quotients, computed by: 


ipixleft = (Xleft[I] — Xmin)/deltaX ; 
ipixright = (Xright[l] — Xmin)/deltaX ; 
We compute the slope of the triangle side: 
slope = (Yright[I] — Yleft[])/CXright[/] — Xleft[l]); 


(Contrary to the previous two assignments, no truncation will take place, since s/ope 
is a floating-point variable.) Proceeding from column ipixleft to column ipixright let 
us assume that we are dealing with column ipix. We are then interested in the row 
number j/ of point I in Fig. 5.17, where the triangle side intersects the boundary 
between the pixel columns ipix and ipix + 1. This row number can be found through 
the screen co-ordinates (ХІ, YI) of intersection point I: 


XI = Xmin + (ipix + 1) * deltaX; 
YI = Yleft(l] + slope * (XI — Xleft[I]); 
jl = (YI — Ymin)/deltaY ; /* implicitly truncated */ 
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Fig. 5.17. Triangle side and pixel columns 
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Let j-old be the old value of jI, computed for column ipix — 1. Then we have 
LOWER{ipix] = min(j-old, jD); 
LOW ipix] = max(j-old, jD; 


where min is defined as the minimum value of its two arguments and max as their 
maximum value. In Fig. 5.17 we һауе LOWER[ipix] = 3, LOW [ipix] = 4. For pixel 
column ipixleft, we assign the truncated quotient (Yleft[I] — Ymin)/deltaY, being the 
row number of the left-hand endpoint, to /.0/4. Similarly, the truncated row number 
(Yright[/] — Ymin)/deltaY of the right-hand endpoint is taken for jZ in the column on 
the extreme right, ipixright. The above method is not yet correct in all cases. If in 
Fig. 5.18 the sides 0 and 1 are dealt with in that order we would obtain the wrong 
result LOWER[ipix] = 6, since the correct value 5, resulting from side 0, would be 
overwritten by the result from side 1. We shall therefore initially assign a very high 
value to the elements of array LOWER and admit these elements only to be 
decremented. It can also be shown that the elements of LOW must be given very 
low initial values and only be incremented afterwards. The values of the arrays 
UPPER and UP are determined analogously. 

For triangle АВС (stored in  TRIANGLE[j] the column range 
ipixmin,...,ipixmax and the array elements LOWER[ipix], UPPER [ipix], 
LOW [іріх], UP[ipix] (ipixmin < ipix < ipixmax) have now been determined, so we 
can now update matrix SCREEN and the linear lists. If the range 


LOW [ipix] + 1, ..., UP[ipix] - 1 


is not empty, and jpix is a row number in this range, pixel (ipix, jpix) is completely 
covered by triangle ABC. We are then interested in point R*, where the line 
through viewpoint E and centre R of the pixel intersects triangle ABC (see Fig. 
5.19). Recall that plane ABC lies at a distance h from E and has n= [a b c]as its 
normal vector. The numbers a, b, с, h are available from TRIANGLE{j]. The 
component of ER in the direction n is EQ. The length EQ is therefore computed as 


— 


ipix Fig. 5.18. LOWER([ipix] depends on two triangles sides 
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Fig. 5.19. Triangle distance 


an inner product: 
ЕО =n- ER 
-[a b c]:[xa yr 1] 
= axg + Був +С 
Since the triangles ERO and ЕК*О* are similar, we have 
ЕК" _ ER 
EQ* EQ 


Hence 
EQ*- ER 
EQ 
hV (xg + ук + 1) 
= axr + Бур +c 


ЕК* = 


The fields tr.cov of all matrix elements SCREEN([ipix |[jpix] are initially set to —1, 
and all fields tr.dist are given the high initial value 10%. For each pixel (ipix, jpix) 
that is completely covered by triangle j these two fields are updated if ER* is less 
than the value already stored in tr.dist. The new triangle distance ER* is then stored 
in field tr.dist and j is written in field tr.cov. 


5.5 AN IMPROVED PROGRAM 


In our program HIDLINPIX most of the new aspects are realized in function 
counter clock. This is a quite complex function. Remember that the main differences 
between our previous program HIDLIN and this new program concern polygons 
and pixels. The decomposition of a polygon into triangles is carried out in essentially 
the same way as in Section 3.5. Every time a complete polygon has been read we 
test whether it is a backface. This test is applied to the triangle whose vertices are 
the first three vertices of the polygon. The test is similar to the one in program 
HIDLIN and is performed in function counter.clock, hence its name. If the first 
triangle turns out to be a backface, the entire polygon is ignored. Otherwise the 
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polygon is not a backface and we have to determine the orientation of three 
successive vertices a great many times. In the same way as in Section 3.5 we thus 
find the next triangle to be cut off the polygon. During this process we want the 
diagonal length too, which is delivered through the fourth argument of the function. 
Finally, the fifth argument is a code to tell the function what we want. The program 
text gives details about this, which we shall not repeat here. 

If a polygon is not a backface not only the triangles but also the edges of the 
polygon have to be stored, since they are the line segments that eventually will be 
drawn, as far as they are visible. Now suppose that two polygons PORS and 
ABCOP occur in the input data and that neither of them is a backface. Since these 
polygons share the edge РО this line segment will be stored twice if we do not 
prevent it. For reasons of efficiency we shall make sure that PQ is stored only once. 
This means that before storing РО we have to search for it in the list of stored line 
segments; should it already be present in this list then we will not store it for the 
second time. АП this has to be done efficiently with respect to time and space, since 
normally there are a considerable number of line segments. There are several ways 
to implement this. Instead of searching, it is attractive to use one of the vertex 
numbers of P and О, say the less of them, as a subscript of an array. As in program 
HIDLIN, we have array VERTEX, where the rectangular co-ordinates of each 
vertex are stored. We now extend each element of this array with a pointer to an 
integer. Actually this pointer will point to the first of a sequence of integers. The C 
functions malloc and realloc enable us to allocate memory space only when we 
actually need it, in a very flexible way. The function realloc turns out to be 
especially useful here. Let us, for example, assume that, in this order, the following 
pairs of vertex numbers are given, each pair denoting a line segment to be stored if 
this has not to be done yet: 


кюбе ©‏ ن رن 
м‏ © © س ت 


This information will be stored as follows: 


i VERTEX(i] 
l x у 2 connect 
1 + d 
0 — 3 2 1 3 
1 — 1 3 
2 — 0 
3 — 0 


The connect field of VERTEX[0], for example, points to the sequence 3, 2, 1, 3. 
The first 3 says that three vertex numbers will follow, namely for the line segments 
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(0, 2), (0, 1) and (0, 3). In C notation we can write 
ptr = VERTEXT[0] : connect; n = *ptr; 


In the example we have п — 3, and 
* (ptr +1) 22 * (ptr +2) =1 * (ptr +3) = 3 


Note that when line segment 3 0 is to be stored this number sequence is converted 
to 0 3 for the sake of uniqueness. If it is not already present, it is then stored in the 
sequence starting in VERTEXT[0] instead of VERTEX[3]. The search for each line 
segment to be stored is now restricted to a usually small sequence, so this method 
will be fast. Since we allocate only memory that we actually need it is also 
economical with respect to space. In the program all this takes place in the function 
add_linesegment. 

To a large extent program HIDLINPIX is similar to program HIDLIN of Section 
5.3. This especially holds for the heart of the program, namely function linesegment. 
In HIDLIN this function had to match the given line segment РО with all triangles. 
їп HIDLINPIX, however, the set of triangles that have to be taken into 
consideration is usually much smaller. Before invoking function /inesegment the 
main program builds this set on the basis of the given line segment РО and the 
linear lists starting in matrix SCREEN. This is done in two steps. First, the arrays 
LOWER and UPPER are used again, but now to denote the pixels that line segment 
РО passes through. This is done in a similar way as for triangles, so we shall not 
discuss this in detail. In the second step the pixels thus selected are used to find the 
triangles that they are associated with. If PO passes through pixel (ipix, jpix), the 
triangles in the linear list starting in SCREEN([ipix][jpix] are added to the set of 
triangles that we are building. The term 'set' is used here to emphasize that the 
array ¢trset will contain each triangle number at most once. Suppose that 
trset[0], ..., trset[ntrset — 1] is the sequence of triangles that we already have. 
(Initially ntrset = 0.) А new triangle number trnr is to be added to the sequence only 
if it does not yet occur in it. This is done by first placing it at the end: 


trset[ntrset] = ntr; 


as a so-called sentinel. We then search the sequence for trnr, starting at trset[0]. 
Since it will certainly be found, we need only test for equality to terminate the loop. 
After the loop is terminated, we test if jnr — ntrset. If so, we know that jtr did not 
yet occur in the sequence before we added it, and we increment ntrset. If jtr < ntrset, 
the triangle did already occur in the sequence, and ntrset is not incremented, which 
means that the sentinel is not effectively added to the set. 

Here is the complete program HIDLINPIX: 


/# HIDLINPIX: A program for hidden-line */ 
/* elimination, using device-independent pixels. */ 


#include <stdio.h> 
#include <math h> 
#include ctype. h> 
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define тах (х, у! CCx )2 (у)? (х): Cyd) 

define тіп(х, у) ((х) (у)? (х): Cyd) 

define тахЗ(х, у, 2) ( (х) > (у) ?тах(х, 2): тах (у, 2)) 
define тміпЗ(х, у, 2) ( (х) < (у) Pmin(x, 2): min(y, 2) ) 
define xwhole(x) CCint) (CCx)-Xmin)/deltax)) 

#define ywhole(y) CCint)CCCu) -Ymin?)/deltaY)) 

#define xreal(i) (Xmint+( (1) )#deltax) 

#define M -1000000.0 


define nvertex 1200 

/* maximum number of vertices #/ 

#define ntriangle 800 

/# maximum number of 

define Nscreen ЗО 

#define big 1. e30; 

define NPOLY 400 

/* maximum number of vertices of a single polygon */ 

define nntrset 200 

/* maximum size of set of triangles associated with 
a single (long) line segment 


(no backface) triangles to be stored */ 


*/ 


int ntr-O, iaux, ipixmin, 
Јріх, top, jbot, j old; 
npoly, isize=sizeof(int), 
LOWCNscreen], 


ipixmax, ipixleft, ipixright, 
1, JI, topcodel3], 
LOWERENscreen]1, 
UPCNscreen], 


ipix, 
POLYCNPOLY], 
UPPERCNscreen], 
trsettnntrsetl, ntrseti 

double vii, vi2, v13, val, v22, v23, v32, 
ерѕ=1е-5, meps=-le-5, oneminus=1-1. e-5 
Xrange, Yrange, Xvp_range, Yvp_range, 
Ymax, deltaX, deltaY, denom, slope 
Xleftt31, Xrightt31, Yleft(31, YrightC3li 


v33, v43, d, ci, 
oneplus=1+1. e-5, 
Xmin, Xmax, Ymin, 


char xmallocí(), #realloc(); 


struct { double 
*pvertex; 


х, Ч, Zi int sconnect; } VERTEX Cnvertex], 


struct € int A, B, Ci double a, b, c, hi } 


TRIANGLE [ntrianglel, #ptrianglei 


struct node < int jtr; struct node next; } #pnodei 
struct < int tr covi double tr disti; 


SCREENCNscreenJ{CNscreen], 


struct node *start; } 
*pointer; 


FILE #fpini 


main(argc, argv) 

€ int i, j: A В, 
*ptr, iconnect, 
trnr, gri 


int argci 
С, Р, @ 
10, 11, 


char #argv[]; 
ii, imin, vertexnr, 
i2, code, count, 


double 
ха, 


xO, 
yA, 
XA, YA, 
х, Y, 
diag, 


yO, 
ZA, 


z0, rho, theta, 
xB, yB, zB, xC, 
XB, YB, XC, ҮС, a, 
хе, уе, ze, 

min diag, Xvp_min, 
fx, fy, Xcentre, Ycentre, Xvp centre, 
xP, ЧР, zP, ха, Ча, 26, ХР, ҮР, XQ, 


X1ft, Xrght, Ylft, Yrght; 


phi, x, 
ЧС, zC, 
b, c, h: 


Ч, 2, 


т, 


Xvp max, Yvp min, Yvp max, 
Yvp centre, 


YG, 
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char chi 


if (argc'=2 ii (fpin-fopen(argvLt11, "r"))-zNULL) 
error("Input file not correctly specified"); 


/* Initialize screen matrix */ 
for (іріх=0; ipix<Nscreeni ipix++) 
for (ypix=Oi jpix&£Nscreen; ур1х++) 
€ pointer-&(SCREENLipixlLjpix1); 
pointer-»5tr cov--1; pointer-Str_dist=bigi 
pointer-»5start-zNULL; 
} 


reflo(&x0); reflo(&y0); reflo(&zO); 

printf("Give spherical coordinates rho, theta, phi of\n")i 
printf("vieupoint E (phi = angle between z-axis and ОЕ) \п"); 
scanf("%1f X1f X1f", rho, &theta, &phidi 

coeff(rho, theta, phi); 

init vieuport(£&Xvp min, &Xvp max, &Yvp min, £Yvp max)i 


/* Initialize vertex array */ 
for (1=0; i£nvertex; i++) VERTEXLil.connectzNULL; 


/* Read vertices */ 
Xmin=Ymin=bigi Xmax=Ymax=-big; 
while (skipbl(), ch=getc(fpin), ch!=’F’ && ch!z-'f') 
X ungetc(ch, fpin); 
reint(£i); reflo(£&x); reflo(&y)i reflo(£z)i 
if (i<O ii i>=nvertex) error("illegal vertex number"); 
vieuing(x-xQ, y-yO, 2-20, &xe, &ye, bie); 
if (ze <= eps) 
X printf("Object point O and vertex Xd lie on ", i)i 
error("different sides of viewpoint E. "); 
} 
X=xe/zei Y=ye/zei 
if (X<Xmin) Xmin=Xi if (X>Xmax) Xmax-X 
if (Y<Ymin) Ymin=Yi if (Y»Ymax) Ymax=Yi 
VERTEXCLil. x=xei VERTEXLil. y=yei VERTEXELi]l. z=ze; 
VERTEXCil. connect = ptr = (int#)malloc(isize); 
if (ptr==NULL) error("memory allocation error 1"); 
*ptr-Oi 


/* Compute screen constants 3/ 
Xrange-Xmax-Xmin; Yrange=Ymax-Ymini 
Xvp range-Xvp max-Xvp min; Yvp range-Yvp max-Yvp min; 
fx-Xvp range/Xrange; fy=Yvp_range/Yrange 
d-z(fxCfy ? fx : fydi 
Хсепёте=0. 5#(Xmin+Xmax)i Үсепёге=0. 5#(Ymin+Yma x); 
Xvp_centre=0. 5#(Xvp_min+Xvp_max)i 
Yvp_centre=0. Ss(Yvp min*Yvp max); 
ci=Xvp_centre-d#Xcentrei c2=Yvp centre-d#Ycentre; 
deltaX-oneplussXrange/Nscreen; 
deltaY=oneplus#Yrange/Nscreeni 
/* Nou we have: Xrange/deltaX < Nscreen #/ 


/* Read object faces and store triangles */ 
while (! isspace(getc(fpin))); 
/* The string "Faces: " has now been skipped #/ 
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while (reint(&i)>O) 
€ POLYCOJ=i;i npoly-i; skipbl(’: 
while (ch=getc(fpin), ch != '&') 
€ ungetc(ch, fpin); reint(&POLYCnpoly++]); 
if (npoly==NPOLY) 
error("too many vertices in one polygon"); 
} 
if (npoly==1) error("only one vertex of polygon"); 
if (npoly==2) 
< add_linesegment(POLYCO], POLYL11): 
continue; 
Y 


if ('counter clock(O, 1, 2, &diag, 0)) continue; 
/* backface */ 


for (i-1; ic-npoly; i++) 
X ii=iZnpoly; code-POLYLiil; vertexnr=abs(code); 
if (VERTEXCvertexnrd. connect-zNULL ) 
€ printf("vertex Xd ", vertexnr); error("undefined"); 
F 
if (codecO) POLYLiil-vertexnr; else 
add linesegment(POLYLi-11, vertexnr); 


} 

/* Division of a polygon into triangles, */ 
/* see Section 3-5: */ 
count=1; 


while (npoly>2) 
€ min_diag=bigi 
for (11=0; і1<про1у; i1++) 
< i0= (ii==0 ? npoly-1 : 11-1); 
i2= (і1==про1у-1 ? О: ii+1); 
if (counter clock(iO, il, i2, &diag,O) 8% 
diag<min_diag) 
< min_diag=diagi imin=ili 
> 
y 
il=imini 
i0= (ii==0 ? npoly-1 : 11-1); 
i2= (ii==npoly-1 ? O : 11+1); 

/* store triangle in array TRIANGLE and in screen lists: */ 
counter clock(iO, ii, i2, &diag, count++); 
npoly--i; 
for (іі=і1; iic-npoly; ii++) POLYLiil-POLYLii-*11i 

+ 
} 
fclose(fpin?)i 


/* Add nearest triangles to screen lists:  */ 
for (іріх=0; ipix<Nscreeni ipix++) 
for (jpix=0; ypix<Nscreeni Jpix++) 
< pointer-&(SCREENLipixlLJpix1); 
if (Gtpointer). tr cov >= О) 
€ pnode-(struct node *)malloc(sizeof(struct пойе)); 
if (pnode==NULL) error("memory allocation error 2"); 
pnode->jtr=pointer->tr cov; 
pnode->next=pointer->start:; 
pointer->start=pnode; 
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/* Drau all line segments as far as they 
for (Р=0; P£nvertex; P++) 
< pvertex=VERTEX+P; /# = &VERTEXECP] #/ 
ptr = pvertex->connect 
if (рт == NULL) continue; 
xP = pvertex->xi ЧР = pvertex-»y; zP 
XP= xP/zP; YP=yP/zPi 


are visible #/ 


= pvertex->zi 


for (iconnect=1i iconnect<=#ptri iconnect++) 


€ Q = #(ptrt+iconnect)i 


pvertex-VERTEX4G; /# = &VERTEXCQ] */ 
xQ = pvertex->xi yQ = pvertex-Syi 28 = pvertex—->zi 


XG7xG/zG; YQ=yQ/7zQi 


/* Using the screen lists, ue shall build the */ 
/* set of triangles that may hide points of PG: #/ 


if (XP<X@ ii (XP--XG && YP<YG)) 


£X1ft-XP; Ylft-YP; Xrght-XG; Yrght= 


YG; + else 


€X1ft-XG; Yl1ft-YG; Xrght=XPi Yrght-YP; } 

ipixleft=xwhole(Xlft); ipixrightzxuhole(Xrght); 

denom-Xrght-X1ft; if (fabs(denom)<=eps) denom=eps 

slope=(Yrght-Ylft)/denomi Jjbot-jtop-yuhole(Ylft)i 

for Cipix=ipixlefti ipix<=ipixrighti ipix++) 

€ if (ipix--ipixright) JjI-yuhole(Yrght); else 

JI7ywhole(Ylft*(xreal(ipix*1)-Xlft)sslope); 

LOWERLipixl2min(jbot, JI); ybot=yTi 
UPPERCipixJ=max(ytop, JI); ytop=Jl; 


ntrset-O; 


for (ipix-ipixleft; ipix<=ipixrighti ipix++) 
for (jpix-LOWERLipix1; yjpix<=UPPERCipixl]i ypix++) 


< pointer=%(SCREENCipixICypixd)i 
pnode= pointer->start; 
while (pnode!=NULL) 
€ trnr- pnode->ytri 


/* trnr will be stored only if it is not yet */ 
/* present in array trset (the triangle set) #/ 
trset[ntrsetj]-trnr; /* sentinel */ 


gtr=0; 

while (trsetCytrl'=trnr) jtr++; 

if (ytr==ntrset) 

€ ntrset**; /* this means that 
if (ntrset==nntrset) 


trnr is stored */ 


error("triangle set overfLOW\n")i 


} 
pnode=pnode->next:; 
} 
} 
/* Nou trsetLOJ, ..., trsetlntrset-1] is the set of */ 


/* triangles that may hide points of PG 


linesegment(xP, ЧР, zP, x@ Ча, 289, 
} 
> 
endgr()i 
} 


0); 


*/ 
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A --------------------------------------------------- #/ 
skipbl() 
€ char chi 
do ch-zgetc(fpin); while (isspace(ch)iicomment(ch)?; 
ungetc(ch, fpin); 
} 
/—--------------------------------------_---—_--—---е---е———- + / 
int comment(ch) char chi 
€ int ki 
if Cch=="(") 
< do k=getc(fpin); while (k != ')' && k != EOF); 


return k==")"; 
} eise return О; 


int reflo(px) double #px; 
€ skipbl()i return fscanf(fpin, "Alf", px); 


int reint(pi) int *pii; 
€ skipbl(); return fscanf(fpin, "ЖЧ", pi); 


add linesegment(P, @) int P, Gi 
€ int iaux, *ptr, іі, ni 
if (Р>@) X iaux-P; P=Qi Q=iauxi } 
/* Nou: P < Q #/ 
ptr=VERTEXCP]. connect; n=*ptri 
for Liisi; ii<=ni ii++) 
if (#(ptr+ii)==@) return; /* G already in list */ 
nti 
VERTEXEP1. connect=ptr=(int #)realloc(ptr, (n+i)#isize); 
if (ptr==NULL) error("memory allocation error 3"); 
X(ptr*n)-2G; #ptr=ni 


int counter clock(iO, il, i2, pdist, code) 
int iO, il, i2, code; double *pdist 


/* code = 0: compute orientation; if counter-clockwise, */ 
/3 compute length of projected diagonal AC */ 
/* code - 1: compute a, b, c, hi store the first triangle */ 


/* code » 1: check if next triangle is coplanar; store it */ 


< int A-abs(POLYLiO1), B=abs(POLYCi1]), C-abs(POLYLi21); 
double xA, yA, zA, xB, ЧВ, zB, xC, ЧС, ZC, т, 
xdist, ydist, zdist, 
ХА, YA, XB. YB, XC. YC, hO, 
DA, DB, DC, D, DAB, DAC, DBC, aux, dist 
xR, yRi 
static double a, b, с, hi 


pvertex=VERTEX+A; 
xA = pvertex->xi yA 


Ш 
il 


pvertex->yi ZA pvertex-2rii 


pvertex=VERTEX+B; 
xB = pvertex-»x; yB = pvertex->yi zB 


pvertex-2rii 
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pvertex=VERTEX+C; 
xC = pvertex-»xi ЧС = pvertex-»yi zC = pvertex->2; 
hO = xA * (yB#zC — yC#zB) — 

xB 3 (yA#zC — yC#zA) + 

XC + (gyA*zB ¬ yB#zA); 


if (code==0) 
if (hOïeps) 
€ xdist=xC-xAi ydist=yC-yAi zdist=zC-zAi 
epdist=xdist#xdisttydist#ydist+zdist#zdisti 
return 1; 
) else return Oi 
/* If hO=0, plane ABC passes through E and hides nothing. */ 


/* If hOCO, triangle ABC is a backface. */ 
/* In both cases ntr is not incremented and the triangles  */ 
/* of the polygon are not stored. */ 


if (code==1) 


< a = yA * (zB-zC) - yB * (zA-zC) + yC * (zA-zBD)i 
b = -(xA + (zB-zC) - xB + (zA-zC) + xC + (zA-1B)); 
с = xA * (yB-yC) - xB * (yA-yC) + xC + (yA-yB); 
т = sqrt(atatbeb+c#c); if (r==0.0) r=eps 
а = a/ri b = b/r; c = c/ri h = hO/ri 


> else if (fabs (a#xC+b#yC+c#zC—-h)>0. 001#fabs(h)) 
€ printf ("Polygon containing vertices Xd Xd Ха ", A, В, С); 
error(" incorrectly specified"); 


} 

if (ntr == ntriangle) error("Too many triangles"); 
ptriangle=TRIANGLE+ntri /# = &TRIANGLECntr] #/ 
ptriangle->A = Ai ptriangle->B = B ; ptriangle->C = С; 
ptrisngle->a = ai ptriangle->b = b i ptriangle->c = с; 
ptriangle->h = hi 


/* The triangle will now be stored in the screen lists #/ 
/* of the associated pixels; first the arrays LOWER, */ 
/* UPPER, LOW, UP are defined: */ 
ХА=хА/:А; ҮА=ЧА/ 2А; 
XB=xB/zBi YB=yB/zBi 
XC=xC/zCi YC=yC/2Ci 
DA=XB#YC-XC#YB; DB=XC#YA-XA#YC; DC=XA#YB-XB#YA; D=DA+DB+DC; 
DAB=DC-M#(XA-XB); DAC-DB-M*(XC-XA)0; DBC-DA-M*t(XB-XC) i 
topcodelOJ=(D#DAB>O); topcodeliJ=(D#DAC>O); 
topcodet21-(D«DBC»0); 
XlefttOJ-XA; YlefttOlJ-YA; XrighttO1-XB; YrighttO21-YBi; 
XleftCiJ=XAi Yleftt121-YA; Xrightt11-XC; YrightC1]=YC; 
Xleft[21-XB; Yleft(21-YB; Xrightt21-XC; Yrightt21-YCi 
for (120; 183; 1++) /# 1 = triangle-side number */ 
if (Xleftt1J»5Xrightt13] ii 
(Xlefttl1--Xrightt11 && Yleft[11»Yrightt11)) 

< aux-zXleftt12J; Xlefttll1-Xrightt11; XrightCl]l=aux:; 

aux-Yleftt11J; Ylefttll1-Yrightt11; YrightCli=auxi 
+ 
ipixmin=xwhole(min3(XA, XB, XC)); 
ipixmax-xuhole(max3(XA, ХВ, ХС) ); 
for (ipix-ipixmin; іріх<=іріхтах; іріх++) 
€ LOWERLipixJ-UPLipix1-2100000; 

UPPERCLipix1-LOMWLipix1--100000; 
+ 
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fo 
< 


} 


/* 
/* 
/* 
/* 


fo 
fo 
+ 


pnode->jtr = ntr; pnode->next = pointer->start; 
pointer->start = pnode 
} 
+ 
ntr++; 
} 
D —— 


HIDDEN-LINE ELIMINATION 


т (120; 1<3; 1++) 
ipixleft-xuhocle(Xleftt11); ipixright-xuhole(Xrightt11); 


denom-Xrightt11-Xleftt11; if (denom==0.0) denom=1. Oe-10; 


slope-z(Yrightt11-Yleftt11)/denomi jy old-yuhole(Yleftt13J 
for (ipix-ipixleft; ipix<=ipixrighti ipix++) 
< if Cipix==ipixright) jI-yuhole(Yrightt11); else 
Jl=ywhole(YleftCli+(xreal(ipix+1)—-XleftCli)#slope)di 
if (topcodeLr11) 
€ UPPERLipix1-2max3(j old, JI, UPPERLipix1); 
UPCipixJ=min3(y_old, JI, UPLipix1); 
} else 
€ LOWERLipixl1-min3(j old, ЈІ, LOWERCipixl); 
LOMLipixJ2max3(J old, JI, LOWLipix1); 


Ji 


*/ 


Y 
J old-JIi; 
} 
For screen column ipix, the triangle is associated only */ 
with pixels in the rows LOWERCipixl,..., UPPERCipix] 
The subrange LOWLipix1*1,...,UPLipix1-1 of these rows 


denote pixels that lie completely whithin the triangle 


т (ipix=ipixmini іріх<=іріхтах; ipix++) 
т (ypix=LOWERCipix]i Jpix<=UPPERCipixli Јріх++) 
pointer=&(SCREENCipixICypixd)i 
if (ypix>LOWLipix] && JjpixcUPLipix1) 
€ xR=Xmin+(ipix+O. 5)#deltax 
yRzYmin-*(jpix-*O. 5)#deltaY; 
denom=a#xR+b#yR+c#di if (fabs(denom)<=eps) denom=eps; 
dist=h#sqrt( xR#xR+yR#yRt1) /denomi 
/* The line from viewpoint E to pixel point (xR, yR, 1) 
/* intersects plane ABC at a distance dist from E. 
if (dist < pointer-»5tr dist) 
€ pointer->tr_cov=ntr; pointer->tr_dist=disti 
} 
} else /* Add triangle to screen list: */ 
{ pnode=(struct node *)malloc(sizeof(struct node)); 
if (pnode==NULL) error("memory allocation error 4"); 


error(str) char #stri < printf("Z%s\n", str); exit(1)i У 


*/ 
*/ 


*/ 
*/ 
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IMPROVED PROGRAM 


eff(rho, theta, phi) double rho, theta, phii 

double th, ph, costh, sinth, cosph, sinph, factor; 

factor=atan(1.0)/45. Oi 

/* Angles in radians: */ 

th-thetas*factor; ph-phisfactori 

costh-cos(th); sinth=sin(th); 

cosph=cos(ph)i sinph=sin(ph)i 

/* Elements of matrix V, see Eq. (4-9): */ 

vil=-sinthi vi2=-cosph#costhi vi13=-sinph#costh; 

v2l=costhi v22=-cosph*sinthi v23=-sinph#sinth; 
v32=sinphi v33--cosphi 

v43-rhoi 


ешіпд (х, yr 2, рхе, руе, рте) 
double x, у, 2, #рхе, #руе, *pze 
/* Ege coordinates, computed as in Eq. (4-2): */ 


*pxe = vil*x + v2al#yi 
*pye = vi2*x + v22*ty + v32*ri 
*pze = у1З#х + vay + v33*z + v43i 


nesegment(xP, ЧР, zP, ха, Ч, 289, КО) 
double xP, ЧР, zP, ха, Ча, zQ; int КО; 
/# Line segment PG is to be drawn, as far as it is not 


/# hidden by the triangles trsetCkO] to trsettntrset-11. 


int у, k-kO, worktodo=1, А, B, C, і, Pbeyond, Qbeyond, 
outside, Poutside, Goutside, eA, eB, eC, sum, imini 
double a, b, c, h, ҺР, hG, ri, r2, r3, 
xA, ЧА, zA, xB, yB, zB, xC, ЧС, zC: 
ЧА, dB, dC, labmin, labmax, lab, mu, 
xmin, ymin, zmin, xmax, ymax, zmax, 
C1, C2, C3, Kl, K2, K3, denomi, denom2, 
Cpos, Ppos, Gpos, aux, epsli 
while (k<ntrset) 
< y=trsetCk]l; 
ptriangle=TRIANGLE+ ji /# = £&TRIANGLEEL JJ s/ 
a=ptriangle->ai b=ptriangle->bi c=ptriangle->c; 
h=ptriangle->h; 


Test 1: */ 

hP=a%xP+b#yP+c#z2Pi hG-a*xG*b*yGecxzG; 
epsl=epsteps#hi 

if (hP-h<=epsl && h@-h<=eps1) {k++;i continue; } 

/* PQ not behind ABC #/ 

Test 2: */ 

Ki=yP#zQ@-yQ#zPi K2=zP#xQ-zQ#xP; K3=xP#yQ—-xQ@#yPi 
A=ptriangle->Ai B=ptriangle->Bi C=ptriangle->Ci 
pvertex=VERTEX+A; 

xA = pvertex-2xi yA = pvertex->y;i ZA = pvertex-»ri 
pvertex=VERTEX+B; 

xB = pvertex->xi yB = pvertex->y; zB 
pvertex=VERTEX+C; 

xC = pvertex->xi yC = pvertex->y; zC 
dA=K1#xA+K2#yYA+KS# ZA; 
dB=K1#xB+K2#yB+K3*2Bi 
dC=K1#xC+K2H#yC+K3#2Ci 


рмет ёех->2; 


pvertex-»zi 


*/ 
*/ 
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/* If ЧА, dB, dC have the same sign, the vertices */ 
/* А, B, C lie at the same side of plane EPG. */ 
eA- dA>eps ? 1 : dA<meps ? -1 : Oi 

eB- dB»eps ? 1 : dB£meps ? -1 : Qi 

eC= dC>eps ? 1 : dC<meps ? -1 : Oj 

sum = еА+еВ+есС; 

if (abs(sum)>=2) { k++; continue; } 


/* If this test succeeds, the (infinite) line РО */ 
/# lies outside pyramid EABC (or the line and the */ 
/* pyramid have at most one point in common. */ 
/* If the test fails, there is a point */ 
/# of intersection. */ 


/* Test 3: */ 
Poutside-Goutside-O; labmin-1.; labmax-O.; 
for (1=0; 1<З; i++) 
< Ci=yA#zB-yB#zAi C2-zA*xB-zB*xA; C3=xA#yB-xB#yAi 
/* C1 x + C2 у + СЗ 1 = О is plane EAB #/ 
Cpos=C1#xC+C2#yC+C3#72Ci 3 
Ppos=Ci#xP+C2#xyP+C3#2P; 
Gpos=Ci#xQ+C2#yQG+C3#7zQG; 
denomi=Gpos-Ppos; 
if (Cpos»eps) 
+ Pbeyond- Ppostmeps; Gbeuyond- Gpos<meps; 
outside- Pbeyond ££& Qpos<=eps i! 
Qbeyond ££& Ppos<=epsi 
+ else if (Сроѕ<терѕ) 
+ Pbeyond= Ppos>epsi Gbeyond- Qpos>eps 
outside= Pbeyond ££& Gpos»-meps |: 
Qbeyond ££& Ppos>=mepsi 
) else outside-1; 
if (outside) breaki 


lab- fabs(denoml)<=eps ? 1.e7 : —-Ppos/denomi; 
/* lab indicates where PG meets plane EAB s/ 
Poutside i= Pbeyond; 

Goutside i= QGbeyond; 

denom2=dB-dA; 

mu= fabs(denom2)<=eps ? 1.e7 : -dA/denom2; 


/# mu tells where AB meets plane EPG #/ 
if (mu>=meps %% mu<=oneplus && 
lab>=meps && lab<=oneplus) 
X if (lab<labmin) labmin=labi 
if (lab>labmax) labmax=labi 
} 
аох=хА; хА=хВ; xB=xCi xC=aux; 
aux=yAi yA-yB; yB=yCi yC=aux; 
aux=zAi zA-iBi zB-zC; zC-auxi 
aux=dAi dA=dBi dB=dCi dC=aux; 
} 
if (outside) {k++i continue; } 


/# Test 4: #/ 
if (!(Poutside ii Goutside)) 
X worktodo=0; break; /* PG invisible */ 
› 
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/* Test 5: */ 
ri-xG-xP; r2-2yG-yP; r3-zG-zP; 
xmin=xP+labmin#rii ymin-yP*labminsr2; zmin=zP+labmin#r3i 
if Ca#xmin+b#ymin+c#zmin-h<-epsi) < k++; continue; ) 
xmax=xP+labmax#ri;i ymax=yP+labmax#r2i zmax=zP+labmax#r3i 
if Ca#xmax+b#ymax+c#zmax-h<-epsi) < k++; continue; У 


/* If this test succeeds, an intersection of PG 


and the pyramid lies in front of plane ABC. */ 
/* Test 6: */ 
if (Poutside ti hPch-epsl1; 
linesegment(xP, ЧР, zP, xmin, ymin, zmin, К+1); 
if (Goutside ii hGch-epsl1) 


linesegment(xG, ЧО, zG, хтах, ymax, тах, К+1); 
worktodo-O; break; 
+ 


if (worktodo) 
< move(d#xP/zP+ci, d#yP/zP+c2) 
drauw(dsxG/zG*ci, d#yQ@/zQtc2)i 


init_viewport(pXMIN, pXMAX, pYMIN, pYMAX) 
double #pXMIN, 3pXMAX, 3pYMIN, 3tpYMAX; 

< double XMIN, XMAX, YMIN, YMAX, len=0. 2; 
printf("Give vieuport boundaries XMIN, XMAX, YMIN, ҮМАХ\п"); 
scanf ("AIF Alf X1f X1f", &XMIN, &XMAX, &YMIN, £YMAX); 
/* Shou the four vieuport corners:  */ 
initgr(); 
move(XMIN, YMIN*1en); drau(XMIN, ҮМІМ); drau(XMIN-*1en, YMIND; 
move(XMAX-1en, ҮМІМ); draw(XMAX, ҮМІМ); drau(XMAX, YMIN+len); 
move( XMAX, ҮМАХ-1еп); drau(XMAX, YMAX); drau(XMAX-1en, YMAX) i 
move(XMIN+len, YMAX); draw(XMIN, YMAX); drau(XMIN, YMAX-1en); 
move ((XMIN+XMAX)/2, YMIN)i draw( (XMIN+XMAX) /2, YMIN) i 
7* Dot in the middle of bottom viewport boundary 

for orientation. */ 

3XpXMINZXMIN; #pXMAX=XMAXi #pYMIN=YMINi #pYMAX=YMAX; 


If we apply this program to the input data at the beginning of Section 5.4, we 
obtain Fig. 5.20 as a result. 
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5.3 


5.4 


HIDDEN-LINE ELIMINATION 


jA 


Fig. 5.20. Letter A in perspective 


EXERCISES 


Compose an input file for HIDLINPIX to draw a simple house in perspective. 
Use HIDLINPIX for one of the Exercises of Chapter 4, but choose a fixed 
value for n, for example n =3. In contrast to Chapter 4, we can now choose 
any viewpoint, without worrying about hidden edges. 

Which modifications of HIDLIN (or HIDLINPIX) are needed if we wish to 
specify the viewpoint in rectangular co-ordinates? 

In this chapter hidden line segments were completely omitted. However, 
engineers often represent them by dashed lines. Investigate how HIDLIN or 
HIDLINPIX can be modified to achieve this. 


СНАРТЕК 6 


Some applications 


6.1 INTRODUCTION 


Program HIDLINPIX is a general program to draw objects in perspective. 
However, despite the good facilities it offers compared with its predecessor 
HIDLIN, a considerable amount of work is still involved in preparing an input file. 
This task can be relieved in two ways. First, there are various types of devices for 
graphic input (see Newman and Sproull, 1979). We shall not discuss these but 
restrict ourselves to device-independent software. Second, we observe that many 
objects have a certain regularity, which enables us to produce the required file 
automatically. Programs for this purpose are executed before the main task, so they 
are sometimes called pre-processors. In this chapter we shall discuss some examples 
of such programs. They all generate a file, which is to be read by the program 
HIDLINPIX. As in Section 5.4, this file has the following structure: 


xO yO zO (Co-ordinates of central object point) 
vertexnumber x y z 


vertexnumber x y z 


Faces: (This keyword is required here) 
vertexnumber ... vertexnumber# (Polygon vertices) 
vertexnumber ... vertexnumber# 


Remarks 


(1) Throughout the file we may insert comment between parentheses (. . .). 

(2) After the keyword Faces the vertex numbers of each polygon are given in 
counter-clockwise order. The final vertex number of each polygon is 
immediately followed by the character #. 

(3) If after Faces a sequence consists of only two vertex numbers, it denotes a 
simple line segment to be drawn (as far as it is visible). This facility enables us to 
draw line segments that are not object edges. 

(4) If a face has a hole we transform it into a polygon by introducing an artificial 
edge. The latter will not be drawn if the second vertex number of the number 
pair for that edge is made negative. For example, if in Fig. 6.1 the inner 
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Fig. 6.1. Hole in rectangle 


rectangle is a hole, we can introduce the artificial edge 2-6 and specify: 
12-65876 -2 3 4# 


Note that if we follow the edges in this order, facing the next vertex, we always 
find the region we are describing on our left-hand side. This means that the 
vertices of the hole are traversed clockwise instead of counter-clockwise. 


In this introduction we shall draw a picture that is a logical extension of the 
example in Section 5.5. Instead of a single solid letter A we shall draw n copies of it 
in a row, where n can be any positive integer. In itself this picture is not likely to be 
of any value, but there is an aspect which we often encounter in practice, namely 
that a rather irregular portion of the picture (a single letter A) has a great many 
duplicates. Then the co-ordinates (and the vertex numbers) of the copies can be 
computed from the original. The following program is a pre-processor. It generates 
the file A. DAT to be read by HIDLINPIX. The program asks for the number of 
letters A to be drawn: 


/* LETTERSA: A Preprocessor for HIDLINPIX #/ 
#include <stdio.h> 
#define thickness 10 
int base; 
FILE Эрт 
main() 
€ int no i; ja пт» xi 
static int Xt201, /* By default initialized to O */ 
yC20] = 4-30, -20, -16 16, 20, 30, О, -12, 12, O}, 
21201 = X О. O, 8, 8, О, O, 60, 16, 16, 40}; 
printf("Hou many letters?\n"); 
scanf("Xd",&n)i 


fp=fopen("a. dat", "u")i 

fprintf (fp, "Af ZF XfNn", —(п-О. S)#thickness, О. О, 30.0); 
for (J=10i )<20; j++) 

€ XCjyJ = -thickness; YCy] = YLj-101; ZCyJ= ZLj-101; 

E 


for (i=Oi itni i++) 

for (J=0; у<20; J++) 

{< nr = 20#i+j; x = —2*is*thickness*XEJ1; 
fprintf(fp, "Xd Xd Xd XdNn", пг, x, ҮСуЈ, ZCyI)i 

} 


INTRODUCTION 133 


fprintf(fp, "Faces: \n")i 

for (і=0; itmi i++) 

{ base = 20%1; 
wi2( O, 1, 2, En 4, 5, 6, -9, 8, 7, 9, bi}; 
wi2(10, 16, -19, 17, 18 19, -16, 15, 14, 13, 12, 11); 
w4( 1, 11, 12, 2); 
w4( 2, 12, 13, 377 
ш4 (14, 4, З. 19) 
ш4( 7, 8, 18, 17); 
w4( 7, 17, 19, 9)i 
w4(18, 8, 9, 19)i 
w4( 5, 15, 16, 6)i 
w4(10, O, 6, 16); 
w4(10, 11, 1, 0); 
w4(14, 15, 5, 4); 

› 

fclose(fp)i 

} 


ш12(а, b, с, d, е, £f, д, hı і, у, К, 1) 
int a, b с, d, е, f. 9, he і, у, К, 1; 
€ fprintf(fp, 
"хаа 44d X4d хаа 74а X4d 74а хаа %4d %4d 74а Х4Ч4#\п", 
s(a), s(b), sc), s(d), sled, 5(#), 
s(g), s(h), 5(1), sj), sk), 5(1)); 
} 


u4(a, b, с, d) int а, b, c, di 
€ fprintf(fp, "%4d %4d 44d %4d#\n", s(a), s(b), sc), s(d)); 
} 


int s(p) int pi 
{ return p>=0 ? p+base : p-base 
} 


To understand this program, see also Section 5.4, where the input file to draw a 
single letter A is given. In contrast to that file, lines of the type 


vertexnumber x y z 


are now written in increasing order of vertex number, but this has no influence on 
the computation. To avoid repeated occurrences of somewhat tedious fprintf calls 
we have introduced the functions w12 and w4. They in turn call the function s, 
which increments the local vertex numbers (—20) by 20i, when dealing with the ith 
letter. Obviously this incrementation is replaced with a decrementation if the vertex 
number was made negative. 

If this program is executed and we type 20 for the number of letters and 700, 80, 
80 for the spherical co-ordinates of E, we obtain a quite extensive file A. DAT. 
Executing HIDLINPIX with this file gives the result shown in Fig. 6.2. The 
computing time on a PRIME 750 was about 17 s. 


(Al 


Fig. 6.2. Output produced by LETTERSA and HIDLINPIX 
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6.2 HOLLOW CYLINDER 


Many objects are bounded by curved surfaces. We can often approximate such 
surfaces by polygons. An example is a hollow cylinder as shown in Fig. 6.3(b). For 
some sufficiently large integer n we choose n equidistant points on the outer circle 
(with radius R) of the top face and n similar points on the bottom face. Then we 
approximate the outer cylinder by a prism, whose vertices are these 2 points. There 
is also an inner cylinder, whose radius is r (r < R). The two cylinders have height Л 
and they share the z-axis of our co-ordinate system as their cylinder axes. The inner 
cylinder is approximated by a prism in the same way as the outer one. The bottom 
face lies in the plane z =0 and the top face in the plane z=h. A vertex of the 
bottom face lies on the x-axis. For given values of n, R, r, h the object to be drawn 
and its position are then completely determined. We shall first deal with the case 
n — 6 and generalize this later for arbitrary n. We number the vertices as shown in 
Fig. 6.4. 

For each vertex i of the top face (1 i « 12) there is a vertical edge that connects 
it with vertex i+12. We can specify the top face by means of the following 
sequence: 


1 2 3 4 5 6 —12 11 10 9 8 7 12 -6# 


Here the pairs (6, ~12) and (12, —6) denote an artificial edge. In Fig. 6.4 the 
bottom face is viewed from the positive z-axis, but in reality only the other side is 
visible. For the bottom face the orientation is therefore the other way round, so that 
we specify: 


18 —24 19 20 21 22 23 24 —18 17 16 15 14 13# 


Since n=6, we have 12 = 2п, 18 = Зп and 24 = 4n, so the above sequences are 
special cases of: 


1... n —2n 2n-1...n+1 2n -n# 


pp 
e 


(a) (b) 
Fig. 6.3. (a) п =6; (D n = 100 
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3 15 
4 9 2 16 21 " 
10 8 22 20 
11 7 23 19 
i 12 À 7 X 13 
6 18 
Top face Bottom face 


Fig. 6.4. Point numbering 


and 
Зп —4n 3n +1 ... 4n —3n 3n —1 ... 2n+1# 


Let us define 


Then the rectangular co-ordinates of the vertices on the top face (vertex numbers 
151,...4 24) are: 

x; = R cos iô 

y=Rsinid (i=1,...,n; outer circle) 

a= h 


x; =rcos(i-n)d 
yJ; =r sin (i -n)ó (i=n+1,..., 2n; inner circle) 
z-Hh 

For the bottom face we have: 


Xj = Xi-2n 
Yi = Yi-2n (i=2n+1,...,4n) 
z,=0 


Here is the program for a hollow prism. By choosing п large enough we obtain a 
good approximation of a cylinder, as Fig. 6.3(b) shows. 


/# HOLLOW CYLINDER (preprocessor for HIDLINPIX) */ 
#include <stdio.h> 
#include <math. h> 
main() 
€ FILE #fp; 

int n, у, К, l, i, mi 

float т, К, pi, alpha, cosa, sina, х, y, Z: 

delta, h, radius, hitei 
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printf("Give number n (n points on a circle): "); 
scanf("Xd", £&n)i 

printf("Give cylinder height: "); 

scanf("Xf", &h); 

printf("Give large radius R and small radius r:"); 
scanf ("AF Xf", SR, т); 


fp=fopen("cyl. dat", "u"); 
pi=4. Otatan(1.0); delta=2. Oxpi/n; 
fprintf(fp, "0.0 0.0 Z7.2f*n", . S#h)i 


for (і=1; i<t=ni i++) 
X alpha=itdeltai cosa-cos(alpha); sina=sin(alpha)i 
for (120; 1<2; 1++) /* 1=0: outer, 1=1: inner circle */ 
< radius- (1==0 ? В: т); 
for (m=O; m<2i m++) /# m=O: top, m=1: bottom boundary */ 
€ kzitlsnems234n; hite= (m==0 ? h : 0); 
fprintf(fp, "Xd 29. SF 49. SF 49. 5f\n", 
kı radius#cosa, radius*sina, hite); 


y 
} 
fprintf(fp, “Faces: \n")i 
/* Top boundary face: */ 
for (i=1i 1<=п; i++) fprintf(fp, "d\n", 1); 
fprintf(fp, "ZdNn", -2s3n)i 
for (i=2#n-1; і>=п+1; i--) fprintf(fp, "Ха\п", i)i 
fprintf(fp, "Ха “~d#\n", 2*n, -n); 


/* Bottom boundary face: */ 

fprintf(fp, "Xd XdNn", Ben, —4*n)i 

for (ї=З#п+1; 1<=4#п; i++) fprintf(fp, "XdNn", 1); 
fprintf(fp, "ХЧ\п", -3#n); 

for (і=3#п-1; 1>=2#п+2; i--) fprintf(fp, "Xd*n", 1); 
fprintf(fp, “"Zd#\n", 2#п+1); 


/* Vertical lines: */ 

for (i=ii 1<=п; i++) 

< fprintf(fp, "Xd Xd Xd “Ad#A\n", 
J=i/n+l, i, i*23*n, у+2+#п); 

fprintf(fp, "Xd Xd Xd Z~Ad#\n", 

itn, jtn: Jj*93*n, i-*3Un)i 

x 

fclose(fp); 


6.3 BEAMS IN A SPIRAL 


Our next example is a spiral as shown in Fig. 6.5. It is built from horizontal beams 
with length /, width w and height w. The bottom beam lies in the xy-plane, as shown 
in Fig. 6.6. Starting at the bottom, each next beam position is obtained by rotating 
the preceding beam about the z-axis through 90° and by incrementing its 
z-co-ordinates by a value w at the same time. We number the beams 0,1,...,n—1 
from bottom to top. Beam i has vertex numbers 8i, 8i+1,..., 8i + 7, assigned in a 
systematic way, such that the vertex numbers of beam 0 are as shown in Fig. 6.6. 
Every point (x', у’, z') of beam i + 1 can be obtained by rotating the corresponding 
point (x, y, 2) of beam i through 90° about the z-axis (and by setting z' =z + w), so 


we have 
cos 90° sin pii 


E у] yl e cos 90? 


137 


BEAMS IN A SPIRAL 


Fig. 6.5. Spiral of beams 


Fig. 6.6. Vertex numbers of beam 0 
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This can simply be written x' — —y and у'=х 
program. 


/# BEAMS (preprocessor for HIDLINPIX) #/ 
#include <stdio.h> 
main() 
+ FILE fpi 
int i, ys п, Ai 
float 1, ш, xA, yA xB, yB, 
printf("Hou many beams?\n"); 
printf("The beam measures 1 


хс, YC: 
scanf("Xd 


. All 


хр, 


SOME APPLICATIONS 


this is used in the following 


yD, a, b; 


ўт); 


aux, Zi 


“ 
D 


X ш x ш. \nGive 1 and ш: "); 


scanf("ZXf Xf", £l, &w)i 
fp=fopen("beams. dat", "u"); 
fprintf(fp, "О.О 0.0 Х#\п", n#w/2.0); 
/# central object point s*/ 
a-O.S9*1; b=a-w; 
хА=а; уА=-а; хВ=а; yB=a; xC=bi yC=ai xD-b; yD=-a; 
for (1=0; 1<п; i++) 
{ for (j20; )<2; J++) 
€ 1=(1+))#ш; A-8*it4*J; 
fprintf(fp, "Xd Xf Af XfNn", А, xA, ЧА, 1)i 
fprintf(fp, "Xd Xf Xf XfNn", A*l, xB; yB, 2); 
fprintf(fp, "Ad AF ZF XfNn", Ate, XC, yC, 2); 
fprintf(fp, "Ad Xf Xf XfNn", А+3, xD, yD, 2); 
+ 
aux=xAji xA=-yAi yA=auxi 
aux-xBi xB--yB; yB=auxi 
aux=xCi xC=-yCi yC=auxi 
aux-xDi xD=-yDi yD=aux; 
+ 
fprintf(fp, “Faces: \n")i 
for (1=0; i<n; i++) 
< А=8%1; 
fprintf(fp,"Xd 4d Жа Ad#\n", А, A*3, А+2, А+1); /*bottom*/ 
fprintf(fp,"Zd Xd Xd “Ad#\n", A*4, A+5, At+G, А+7); /#top #/ 
fprintf(fp,"Xd Xd Xd Ха#\п", А, A*1, A+5, А+4); /*front */ 
fprintf(fp,"Xd Xd Xd Ха#\п", A*3, A*7, A+6, &A*2);/*back #/ 
fprintf(fp,"Xd Xd Xd Ad#\n", A, A*4, A*7, A+3); /*left #/ 
fprintf(fp,"Xd Xd Xd Xd&*Nn", A*1, A*2, A*6, A+5)i /*right */ 
+ 


fclose(fp)i 


6.4 WINDING STAIRCASE 


The idea of beams in a spiral leads to our next example, a winding staircase, as 
shown in Fig. 6.7. If we mount the stairs on each step our height increases by Л. 


There are n stairs, so the total height is H = nh. 


Each stair has eight vertices, locally 


numbered 0,...,7. We add the endpoints 8 and 9 of a vertical bar to them. This 
bar serves to attach a hand-rail to the stairs. The bars and the hand-rail are very 
thin, so we draw them as line segments. The centre of the staircase is a cylindrical 
pole with diameter 2r. Its axis coincides with the z-axis of our co-ordinate system. 
The hand-rail and the vertical bars are at a distance R from the z-axis (R >r). The 
stairs connect the pole with the vertical bars. The bottom stair is shown in Fig. 6.8. 


Each stair is a beam with length R — r, width 1. 


5h and height 0.2/4. Together the л 


Fig. 6.7. Winding staircase 


Fig. 6.8. Bottom stair 
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stairs constitute a full revolution, so the angle of rotation for a single step is 


m: 
n 

From bottom to top we assign the numbers 0,1,...,n—1 to the stairs. Stair i 
can then be obtained by rotating stair 0 about the z-axis through the angle œ = iô 
and by raising it to the height ih at the same time. Thus each vertex (x, y, z) of stair 
i can be computed from the corresponding vertex (X, Y, Z) of stair 0 as follows: 


cosa sin si 


b у]=[Х vil 
z=Z+ih 


—sina COS а 


In the program that follows the vertex co-ordinates of the bottom stair are stored 
in the arrays X, Y, Z. For the stairs, including the vertical bars and the hand-rail, 
we use the vertex numbers 0,..., 10" — 1. Setting M = 10n, we assign the next n 
numbers M,..., M +n — 1 to equidistant points on the bottom circle of the central 
pole. Finally the integers M +n,..., M +2п— 1 are assigned similarly to the top 
boundary of the pole. 


/* MINDING STAIRCASE (preprocessor for HIDLINPIX) */ 
include zstdio h> 
#include ‘math h> 
main() 
{ FILE *£fpi 
int i, je Mr Ke M, li 
float г, R, pi, alpha, cosa, sina, х, ц, 2, 
delta, һ, Н, XC10], Yt101. ZEL101; 
printf("Give number n (to draw n stairs)\n"); 
scanf("Zd", £n) 
printf("Give height of a single stair step "); 
scanf( "ZF", £h»i; 
printf("Give large radius R and small radius r\n“), 
scanf ("AF Xf", YR, rhi 
fp-fopen("uinding.dat", "ш"), 
p1=4. O#atan(1.0); delta=2 O#pi/n; 
fprintf(fp, "O.O 0.0 Xf*n", О Ssnsth) 
for (Ј=0; J4; ++) ZE JJ=0; 
for (J=4i 58; J++) ZLj1-h/5 
XCOJ=XC1I=XC4I=XCSI=R, 
XC2)=XC€3I=XC6I=XC7I=ri 
YCOI=YC4I3=YCO3I=YC7I=—-. 75*h; 
ҮС1Ј=ҮГ5Ј=ҮГ2]=Ү[Г6]Ј=. 75#h, 
XCBI=Ri YCBI= О, ZCBI=h/10; 
XC9J=R; Yt91-2 О; ZI(91-5*h; 
M=10#n; H=n#h, 
for (1=0, itn, i++) 
t alpha=i#delta, cosa=cos(aipha); sina=sin(alpha)i 
for (у)у=0; Jjz10; J++) 
€ k=10#1+y, 
x=XCjJl#cosa-Yljl#sina; 
y=XCjl#sina+Yljl#cosa; 
z=Z0j14+1 8h; 
fprremtf€fp, “Ad Ze Wf AFN", Ks X2 Ge 299 


TORUS 


х=г#соѕа; y=r#sina; 
fprintf(fp, "Ad Xf XP ARAN, Mei, x» y: O о), 
fprintf(fp, "Xd Xf Xf Ж#\п", М+п+1, x» у, H+5#h); 

> 

fprintfifp, "Faces: Nn"); 

for (1=0i ini i++) 

€ к=10%1; 
fprintf(fp, "Ad Ad Xd Ad#\n", К, k+l, k+5, К+4); 
fprintf(?p, "Ad Xd Xd AdHAN", k*2, k+3, k+7, kt&)i 
fprintf(fp, "Ad 4d Ad Ad#\n", k+l, k+2, k+6, k+5); 
fprintf(fp, "Ad Ad 4d Ad#\n", k*3, k*O, k*4, К+7); 
fprintf(fp, "Zd Zd 4d XdWNn", k*7, k+l, k+5, k+6): 
fprintf(fp, "Xd Xd *d XdfNn", k+4, k*5, k+6, k+7); 
fprintf(£p, “Ad 4d Xd Xd&Nn", k*i. k*O, k+3, k+2):; 
fprintf(fp, "Xd XdÉiNn", k+B, k*9); 
if (izn-1) fprint£(fp, "Ad Ad#\n", k+9, k+19); 

} 


for (і=0; ini i++) 
€ fprintf(£p, "Ad Ad Ad XdWHNn", M+i, M+(i+1)%n, 
M*n-(Citi)Zn, M+n+i); 


} 

for (120; 122; 1++) 

{ for (і=0; ini i++) fprintf(fp, " Xd", M+l#n+i); 
fprint^(fp, "#\п"); 

} 

fclose(fp); 


6.5 TORUS 
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In the winding staircase some dimensions (or rather proportions) were chosen 
arbitrarily. We shall now discuss some examples where all vertices are computed 


from a very limited amount of data. The first is a torus, as shown in Fig. 6.9. 


Fig. 6.9. Torus 
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Fig. 6.10. Basic circles of torus 


The input of the program will consist of three numbers, n, R, r (R >r). In Fig. 
6.10 the large horizontal circle is the centre circle of the torus; the radius of this 
circle is R. We choose 7 equidistant points on this circle as centres of small vertical 
circles with radius r. A parametric representation of the large circle is 


x = R cos « 
y-Rsin« 
z=0 


The point corresponding to а = 0 is the centre of the small circle 


х= К + т со В 
у= 0 


2 —rsinf 


which is also shown in Fig. 6.10. By rotating this particular small circle about the 
z-axis through angles а = ід, where i—1,...,n —1 and ó =2x/n, we obtain the 
remaining n — 1 small circles. On the basic small circle in Fig. 6.10 we choose n 
points and assign the vertex numbers 0, 1,...,n — 1 to them: the point obtained by 
choosing parameter Û = jô is given vertex number j (j = 0, 1, ..., n — 1). The next n 
vertices, numbered n,n+1,...,2n—1, lie on the neighbouring small circle, 
corresponding to i = 1, and so on. In general, we have the vertex numbers i: n + j 
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((20,1,...,n —1,]250, Gases = 1). A rotation through the angle а =i about 
the z-axis is written 


m vil cos а abel 


—sin а COS а 


In our situation the basic small circle lies in the xz-plane, so y = 0, which reduces 
this matrix product to 


X’ =x cos a 
y' = x sin а 


This result could also have been derived immediately from Fig. 6.10. The 
following program generates the file for the torus. 


/* TORUS (preprocessor for HIDLINPIX) */ 
#include <stdio.h> 
#include <math. h> 
main() 
€ FILE #fp; 
int i, yo mi 
float т, Н, pi, alpha, beta, cosa, sina, 
х, х1, yl, z1, deltai 
printf("Give number n (to draw ап п x n torus)\n"); 
scanf("Xd", ®п); 
printf("Give large radius К and small radius т\п"); 
scanf ("AF Xf", YR, Er); 
fp-fopen("torus.dat", "u"); 
piz4.O0*atan(1.0); delta=2. Oxpi/ni 
fprintf(fp, "О. О 0.0 О. О\п"); /# central object point #/ 
for (i=Oi i<ni i++) 
X alpha=itdeltai cosa=cos(alpha)i sina=sin(alpha)i 
for (y=Oi j&ni J++) 


< beta=y#deltai x=R+r*¥cos(beta)i /*y = O */ 
xi=cosa#x; yi-sina*x; zl-r*sin(beta); /* 21 = z */ 
fprintf(fp, "Xd Xf Xf KEN", itn+y, xl, yl, 210); 

+ 


y 
fprintf(fp, "Faces: n"); 
for (1=0; 1<п; i++) 
for (Jj20; JEni J++) 
€ fprintf(fp, "Xd Xd Ad Zd#\n", 
i*ntj, (Citl)ZAn*n*tj, (Citi)Xnnt(CJ*ti)An, isnt(CJtl1)0Zn)i 
+ 
fclose(fp); 


6.6 SEMI-SPHERE 


Figure 6.11 shows the lower half of a sphere. We let the origin of the co-ordinate 
system coincide with the centre of the sphere. Since the whole picture will 
automatically be as large as the viewport, the absolute size is irrelevant, so we can 
choose a radius of unit length. Thus all points (x, y, z) of the semi-sphere satisfy: 


S 


—1<z<0 
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D وج‎ 


Fig. 6.11. Semi-sphere 


Point (0, 0, —1) is given vertex number 0. We shall divide a right-angle into n equal 
angles д; this number п is the only number that our program will read. The relevant 
points on the semi-sphere with x > 0 and y =0 are numbered 1,2, .. . , n, counting 
from bottom to top. Their neighbours with y >0 are numbered n +1, n + 
2,...,2n, and so on. This way of numbering is illustrated by the following table, 
where each (horizontal) row corresponds to points on a horizontal circle; similarly, 
each column corresponds to points on a quarter of a vertical circle: 


i > 0 1 2 4n —1 
J 
шаа СЕ 
| 
п п 2п 3n 4n? 
2 2 n+2 2n +2 (4n — 1)n +2 
1 1 n+1 2n +1 (4n — 1)n +1 


The lowest layer of the semi-sphere consists of triangles, all having point 0 as a 
vertex. The remaining п — 1 layers consist of quadrangles ABCD; each quadrangle 
has two edges AB and DC that are horizontal and mutually parallel. Note that this 
semi-sphere differs fundamentally from the other examples that we have seen. It is 
not a solid object but a surface, which is visible from either side, depending on the 
viewpoint. Since each quadrangle has two faces, we specify both sequences ABCD 
and DCBA in the following program. 


FUNCTIONS OF TWO VARIABLES 145 


/* SEMI SPHERE (preprocessor for HIDLINPIX) */ 
include <stdio.h> 
#include <math. h> 
main() 
€ FILE *fpi 
int i, у, m А, В, С, D, P, Gi 
float pi, alpha, beta, delta, cosa, sina, cosb, sinb; 
printf("Give number n\n"); scanf("Zd", &n); 
fp=fopen("semi. dat", “w“); 
рі=4. Otatan(1.0); delta-pi/(2*n); /* n * delta = pi/2 */ 
fprintf (fp, "0.0 0.0 -O. 5Nn"); /* central object point #/ 
/* R = 1i sphere centre in О #/ 
fprintf(fp,"O 0.0 0.0 1. О\п"); /* first point #/ 
for (1=0; 1<4#п; i++) 
< alpha-i*delta; cosa=cos(alpha)i sina-sin(alpha); 
for ()=1; у<=п; j++) 
< beta-j*delta; cosb=cos(beta)i sinb=sin(beta)i 
fprintf(fp, "Ха Xf Xf Xf*n", 
n#i+y, sinb*cosa, sinb*sina, -cosb)i 
+ 
У 
fprintf(fp, “Faces: \n")i 
for (1=0; i<4#ni i++) 
€ Р=1ї#п+1; @=(1+1)7Ж(4#п)#п+1; 
fprintf (fp, "Ad %d “~Ad#\n", О, P, Q); 
fprintf(fp, "Ха Xd “~Ad#\n", О, G, Р); 
for (J=li yeni J++) 
< A=P+j-1; B=Q+j-1i C=Q+j; D=P+ jJi 
fprintf(fp, "4d Xd Ad Z~Ad#\n", А, В, C, р); 
fprintf(fp, "29 Xd Xd “Ad#\n", D, C, B, А); 
> 
} 
fclose(fp); 


6.7 FUNCTIONS OF TWO VARIABLES 


Program HIDLINPIX was designed to produce perspective drawings of solid 
objects. In two respects we have already seen that we can use it for other purposes 
as well. First, ‘loose’ line segments such as portions of co-ordinate axes can be 
drawn if we include their endpoints in the input file under the keyword ‘Faces’. 
Second, in Section 6.6 we dealt with the surface of a semi-sphere as another 
mathematical abstraction. We shall now go a step further in using HIDLINPIX for 


purposes beyond our original goal. 
We sometimes require a graphical representation of functions of two variables: 


z — f(x, y) (6.1) 


We wish to specify a rectangular domain: 
X min sxs X max 


Ymin = y = Ymax 


In principle, function f can be any continuous function of two variables. As an 


example we shall use the quadratic function 


f(x, у) = ах + by + cxy + dx c ey +g (6.2) 
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(We do not use the letter f as a coefficient because f is the function name.) The 
program will ask for the coefficients a, b, c, d, e, g and also for the integers N, and 
N,. These two integers are used to compute the following length and width of 
elementary rectangles: 


= Xmax — Xmin 


A 
x N. 


Ymax — Ymin 
LN и 


y 


The corners of these rectangles are the grid points (x, y) for which we actually 
compute z=f(x, y). The three-dimensional points (x, y, 2) thus obtained are 
connected by straight line segments parallel to either the xz- or the yz-planes. It is 
these line segments that will be drawn. Figure 6.12 shows the function 


f(x, y) =0.1x? — 0.4у? 


where we have chosen 


—5*x*5, —2<y <2, N, = 20, N,=8 
р = 20, Ө = 50 ф = 80 
We associate a pair of integers (i,j) with each grid point (= 0,..., M; j= 
0, ... , Ny). The point (x, y, z) of the surfaces that corresponds to grid point (i, j) is 


assigned vertex number k =]. (N, +1) +i +1. Then we have: 
= es FIA 
Y = Ymin Fj ` Ay 
z = f(x, у) 


The vertex numbering is illustrated in Fig. 6.13(a), where N, = 3 and Ny = 2. 

In Fig. 6.13(a) we view the surface from the positive z-axis. For example, points 
1, 2, 5, 6 are points on the surface. Unfortunately, these four points will in general 
not lie in the same plane, so we cannot use them as vertices of a polygon. If we 


Fig. 6.12. A quadratic function of two variables 
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(a) (b) 


Fig. 6.13. (a) Points on surface; (b) two triangles 


connect point 1 with point 6 we have two triangles that solve our problem. For the 
general case this is shown in Fig. 6.13(b), where /=k+N,+1. Though not 
co-planar, these two triangles can be used as the required polygons, provided that 
we prevent edge (К, / + 1) from being drawn. Depending on our viewpoint, either 
side of the triangles can be visible, so under Faces we have to specify the triangles 
twice: 


k —(L-1) К+1# (lower-right, clockwise) 

k+1 [+1 —k# (lower-right, counter-clockwise) 
k —(1+1) l# (upper-left, counter-clockwise) 
l [+1 —k# (upper-left, clockwise) 


The minus signs will prevent line segments (К, / + 1) from being drawn. 

Besides the function surface we shall also draw portions of the positive 
co-ordinate axes, as far as they are visible. Their lengths will be specified by the 
user. Finally, the program will also ask for some central z-value to be able to 
determine the ‘central object point’ required by HIDLINPLX. This is not particu- 
larly critical, so a rough estimation will do. 

For quadratic functions of two variables the following program is quite general. 
For other functions f(x, y) we can replace the function f at the bottom of the 
program and remove the program elements that have to do with the coefficients a, 
b, 6; d'ê 2 


/* FUNC: Perspective plot of a quadratic function #/ 
#include <stdio.h> 
float а, b, с, ds е, gi 
main() 
{ FILE fpi 
int i, Je Nx, Ny; ke 1; 
float xmin, xmax, ymin, ymax, hx, һу, x, y, FO), 
zc, хахіѕ, yaxis, zaxisi 
fp=fopen("FUNC. DAT", "ш"); 
printf("f(x, y) = а. х.х + boy. y + c. x.y + а. х + e.y + gin"); 
printf("a, b, c, d, e. 9: "di 
scanf (“AF Xf Xf Xf Xf Xf", Ba, Bb, Bec, &d, Be, &g)i 
printf("xmin, xmax, ymin, ymax:"); 
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scanf ("Xf 
printf("Nx, 


scanf("Xd Xd", 


hx= 


Cxmax-xmin)/Nx; 


printf("Central z-value:")i 
printf("Length of positive axes to be drawn (x, y, 


scanf("Xf", 


Xf Xf Xf", &xmin, &xmax, &ymin, &ymax); 
Ny: ")i 
&Nx, &Nuy); 
hy= (ymax-ymin)/Ny; 


&zc)i 
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27: "M 


scanf("Xf Xf Xf", &xaxis, &yaxis, &zaxis)i 

х= (хтміп+хтах) /2; y= (утіп+утах ) /2; 

fprintf(fp, "Xf ZF XfNn", x, y, ze; 

for (1=0; i<=Nxi i++) 

for ()=0; J<=Nyi j++) 

< xzxmintishx; y=ymin+Jy#hyi 

. fprintf(fp, "Ad Xf ZF XfNn", Jg*ONx*ti)titi, x, y, #(х,у)); 

r 

kz2(ONx*1)3 (Ny *1); 

fprintf(fp, "Xd Xf Xf XfNn", ++k, O., O., O. ); 

fprintf(fp, "Xd Xf Xf ZXfNn", ++k, xaxis, O., O. ); 

fprintf(fp, "Ха Xf Xf XfNn", ++k, O., yaxis, О.); 

fprintf(fp, "Ха Xf Xf XfNn", ++k, O., O., zaxis); 

fprintf(fp, "Faces: n"); 

for (і=0; i<Nxi i++) 

for ()=0; y<Nyi J++) 

€ к=)#(Мх+1)+1+1; 1=k+Nx+1; 
fprintf(fp, "Xd Xd Ad#\n", К, —(1+1), К+1); 
fprintf(fp, "Xd Xd ЖЧ#\п", k+1, 1+1, —К); 
fprintf(fp, "Xd Xd Xd$W Nn", К, —(1*1), 1); 
fprintf(fp, “Ad Xd XdWNn", 1, 1+1, —К); 

+ 

к=(Мх+1)#(Му+1); 

fprintf(fp, "Xd ХЧ#\п", k+1, k*2)0; /+% x-axis */ 

fprintf(fp, "Xd Ad#\n", k+1, k*3); /* y-axis */ 

fprintf(fp, "Xd  Xd&Nn", k+1, kt4);  /* z-axis */ 


fclose(fp); 


У 


float f(x,y) float х,у; 
€ return a*xsx*bsystytcitxitytdstxtestytgi 


› 


6.1 


6.2 
6.3 


6.4 


6.5 


EXERCISES 


Write a pre-processor to draw several pyramids, some of which partly hide 
others, depending on the viewpoint. 

Write a pre-processor to draw a full sphere. 

Write a pre-processor to draw several semi-spheres, some of which partly hide 
others. 

Choose one of the Exercises of Chapter 4 and write a pre-processor for it. We 
are now in a position to give a solution that is general in two respects: 


— Hidden lines are eliminated by computation (in contrast to what we did in 
Chapter 4). 
— The number n is variable. (In Exercise 5.2 we chose n = 3.) 


Write a general rotation program in the sense of Exercise 3.2. The user will 
specify a vector AB and an angle a. The program is to read an input file for 
HIDLINPIX and to write another. In the latter file the co-ordinates will differ 


EXERCISES 149 


Fig. 6.14. A cube of cubes 


from those in the former, according to the given rotation, so that a picture of 
the rotated object will be the result. 

6.6 Write a pre-processor to draw a great many cubes that are placed beside, 
behind and above each other (see Fig. 6.14). 


APPENDIX 


A brief introduction to the C language 


This Appendix does not deal with the complete C language. It is intended for those 
readers who are familiar with another modern programming language and who want 
a brief survey of the C elements used in this book. Anyone who is going to use C 
actively is strongly advised to consult a textbook such as The C Programming 
Language, by Kernighan and Ritchie (1978). 


A.1 BASIC DATA TYPES 


Before we can use any variables we have to specify their types, or, in technical 
terms, we have to declare them, as, for example, in: 


float xA, yA, p, q, т = 5.0, epsilon = 1e — 6; double dd; 
int i; short int si; long int li; 
char ch =° A; 


In their declarations, variables can be given initial values, as is done here for the 
variables r, epsilon and ch. We say that they are initialized. The numbers of bits that 
are used for the various data types depend on the hardware. On many machines, 
type float (single precision floating point) can effectively be replaced with type 
double (double precision floating point) to achieve greater precision. In the 
following we shall not always mention double explicitly, but regard it as a special 
case of float. Conversions from float to double and vice versa will never have effects 
different from what we expect. There is no ‘Boolean’ or ‘logical’ type in the C 
language. The integer value 0 acts as false, the value 1 (or any other non-zero value) 
as true. 


A.2 SOME STATEMENTS 
Consider the following fragment of a program: 


if (x> —0.5) (i = 10; j = 20;) else {k = 30; l = 40;} 

if (u<3.0) {v = 1.8; w =3.4;} 

if (a<b)m=n= 100; 
The meaning of these three if-statements will intuitively be clear. Notice that if is 
always followed by an expression between parentheses. An else-part is optional. 
Assignment statements such as 

i — 10; 
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(and all other statements too) can be grouped together into a single compound 
statement by braces. Thus 


(i = 10; j = 20;) 
is a compound statement. We need them in if-statements if the execution of more 


than one statement has to be dependent on the condition. For example, if we had 
written 


if (u <3.0) v = 1.8; w = 3.4; 
then 
и = 3.4; 


would have been executed anyhow, since the latter assignment statement has 
nothing to do with the condition u < 3.0. Braces are not used if only one statement 
is to depend on the condition, as the last of the three given if-statements shows. 
Here the multiple assignment statement 


m — n — 10; 


assigns the value 10 to both variables m and n. 
A primitive loop can be set up by means of an if-statement and a goto-statement, 
as in 
s=0; i-l; 
again: if (i< =n) {s+ =i; i +; goto again;} 
Here s + = and i + + are short-hand for = +i and i =i + 1, respectively. These 
two lines are not meant as an example of well-structured programming. They merely 
serve to explain two new language concepts. The first is the while-statement, used 
in: 
5 —0; p=]; 
while (i<=n){s+=i;i++;} 


This loop has the same effect as the previous one, so after i is incremented it is again 
compared with п and so on. 
Another important construct is the for-statement, used in: 

$0; 

for (i=l;i<=n;i++)s+=i; 
This piece of program has the same effect as the two previous ones, so in each of the 
three cases the following value is computed. 

s=1+2+...+n 


Should лп have the value 0 (or less than 0), then s will also be 0, since in all three 
cases the test for loop termination is executed at the beginning of the loop. This is 
not the case in the do-while-statement, appearing in: 

i=l; 5 = 0; 

do {s+=i;i++;} while (i< = n); 


Here the test for loop termination is executed at the end, as the notation suggests. 
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Thus the value 1+2+...+n is given to s only if п is positive. Should n be 0 (or 
less), then s obtains the value 1, in other words the inner part of this loop is 
performed at least once. 
The statement 
break; 


can be used to terminate a while-statement, a for-statement or a do-while-statement 
unconditionally whenever this is necessary. Thus the effect of 


i-1; 
while (i< =n) (s + =i; if (s > MAX) goto ready; i+ +;} 
ready: 


is obtained without a goto-statement by writing 


і= 1; 
while (i< =n) {s + =i; if (5 > MAX) break; i+ +3} 
Or 
for (i=1;i<=n; i++) {s+=i; if (s > MAX) break; 


To go immediately to the test for loop termination we can use the statement 
continue ; 
For example, the while-statement 
while (a < b) {a + +; if (a<p){a+=2;b+=1;}} 
can be replaced with 


while (a < b) {a+ +; if (а> = p) continue; a+=2;b+=1;} 


A.3 OPERATORS AND EXPRESSIONS 


Characters such as +, —, < are used as operators to build expressions such as 
a+b-—c<d 
If we write 
X-ty*z 


we require the product y*z to be added to x. Indeed, the multiplication operator * 
has a higher priority than the adding operator +. Operators that have the same 
priority usually associate from left to right. For example, the expression 


a—b+c-d 
means 
(а=) +) а 


There are also operators that associate from right to left. An example is the 
assignment operator —, since the assignment statement 
i=j=k=30; 
has the same effect as 
i = (j = (К = 30)); 
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Note that — must not be used to test whether two values are equal; for the latter 
purpose we have the operator = =. 

We shall now list all C operators; the meaning of some of them will be clear 
immediately, others will be discussed later. The operators are listed in order of 
decreasing priority, but those between a horizontal line and the next one have the 
same priority; they also associate in the same way. 


( ) Function call left-to-right 
[ ] Array element 

Structure member 
-> Structure member using pointer 


(The following operators are said to be unary: they have only one operand. It 
should be noted that the characters —, &, * are also used as binary operators.) 


! Logical not (see below) RIGHT-TO-LEFT 
~ One’s complement 

= Negative 

++ Increment 

== Decrement 

& Address 

* Indirection 

(type) Cast (see below) 

sizeof Size in bytes 

* Multiply left-to-right 
/ Divide 

% Remainder 

+ Ааа left-to-right 
= Subtract 

« Left shift (see below) left-to-right 
>> Right shift (see below) 

« Less than left-to-right 
<= Less than or equal 

> Greater than 

>= Greater than or equal 

== Equal left-to-right 
= Not equal 

& Bitwise and left-to-right 


^ Bitwise exclusive or left-to-right 
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| Bitwise or 


(see below) 
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left-to-right 


&& Logical and (see below) left-to-right 
| Logical or (see below) left-to-right 
di Conditional (see below) RIGHT-TO-LEFT 
= = f= کوټ‎ += == &= ъ= &= A= |=RIGHT-TO-LEFT 
Assignment 
$ Сотта (see below) left-to-right 


An expression may not only yield a value but also perform an action that changes 
the environment. The following line, for example, is an expression: 


i-icl 


This is an assignment expression, which increases i by 1. Less obviously, this 
expression also has a value, namely the new value of i. Hence it makes sense to 
write 


j=S*(i=i+1)+2; 
which, incidentally, can be replaced with 
j=5*(++i)+2; 
Both + +i and i++ have the effect that i is increased by one. However, they 


differ in the value they yield in the following sense. 


+ + i means: increment i and use its (new) value; 
i + + means: use the (old) value of i, then increment 
this variable 


Thus after the execution of 


we have i 2j =6, m=6, п = 5. 
An assignment expression and a semicolon (;), in that order, form an assignment 
statement. Thus 


=ї+ 1; 


is not an expression but a statement (and so is i + +;). If several actions are required 
in a context where only an expression (and no statement) is allowed we can use the 
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comma operator. Consider, for example, the loop 
while (i 2j,j-—,k-i*2*j, k >0) j*=2; 
This works as follows: 
again: i+ =j; j——;k=i+2*j; if (k >0) (j*—2; goto again;} 
(the meaning of j*= 2 is j = j*2, j – = 3 means j =] — 3, etc.) 


The pair ?: forms another unconventional but useful operator. In the so-called 
conditional expression 


cond ? expr1:expr2 


first the expression cond, denoting a condition, is evaluated. Then either exprl or 
expr2 is evaluated, depending on whether the value of cond is non-zero or zero, 
respectively. Remember that true is 1 and false is 0 in the C language. Thus instead 
of 


if (x > y) max = x; else max = y; 


we can write 
max = (x »y?x:yy 


It is important to distinguish between logical and bitwise operators. Logical 
operators are often used, as, for example, in 


(x <0|| y 20&&y «1)&&!(z <0) 
which should be read as 


((x <0) or ((y > 0) and (y < 1))) and not (z <0) 


For the logical operators && and || it is guaranteed that the second operand is 
evaluated only if the first is not sufficient to decide what the final result will be. 
There is therefore no danger that 


i»0&&jli ^k 
should cause a division by zero. 
The result of such a logical expression is always 0 or 1. On the other hand, if we 


use bitwise operators we think in bits, although we use integer variables, as, for 
example, in: 


i5; jS k=i|% 


Here the operations ‘left shift’ (<<) and ‘bitwise or’ (|) are used. Written in binary, 
the values of i, j, К are now as follows. 


0...000101 (= 5) 
0...010100 (=20) 
0...010101 (= 21) 


Il 


i 
j 
k 
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The cast operation can be used to force type conversion. Thus the value of 
(int) 3.95 


is 3 of type integer. 

If both i and j are of type integer, i/j is the truncated quotient. Thus 8/3 has the 
value 2, of type integer, even in a context where float is expected, as in x — 8/3, 
where x has type float. If at least one of the operands a and b has type float, the 
quotient a/b has type float and is not truncated. Float constants contain a period or 
the letter e (or E); integer constants have not. Thus: 


1.5e3/50 has the value 30.0 
7/4.0 has the value 1.75 


7/4 has the value 1 
The remainder operator % can be used for integer operands only: 
37%5 has the value 2 
37% 5.0 is illegal 


(int) 37.9% 5 has the value 2 


In the last example the cast operator (int) has priority over the remainder 
operator %. 

If a float value is assigned to an integer, truncation will take place. Thus the 
statement 


i=3.9; 


assigns the value 3 to the variable i. Even if variable x has type float, the following 
statement will assign the truncated value 7 (converted to float 7.0) to x: 


x = 39/5; 


A.4 LEXICAL TOPICS AND PROGRAM STRUCTURE 


In most cases blank spaces and transitions to new lines have no effect on the 
meaning of the program. Sequences of the form 


/*... */ 


are ignored by the compiler; they are pieces of comment for the human reader. 
Identifiers such as names of variables consist of letters and digits; the first 
character must be a letter. Here the underscore (_) is regarded as a letter. Taking 
this into account, and distinguishing capital and small letters, we count 53 distinct 
letters. 
We can also use an identifier for a constant, as, for example, in 


#define MAXIMUM 1000 


After this ‘pre-processor control line’ we can use the identifier MAXIMUM as just 
another notation for 1000. It is the first example of a macro. A more interesting one 
is MAX, defined by 


#define MAX(x, y) x » y?x:y 
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This line has the effect that, throughout the rest of the program, any string 
MAX (a, b) 


is automatically replaced with: 
a>b?a:b 


To obtain correct results in more complicated cases, it is wise to use parentheses, 
so we should actually define this macro as follows: 


#define MAX(x, y) ((х) >(y)?(x):()) 
There are also control lines for file inclusion. We often use 
#include (stdio.h) 


The effect is the replacement of this line with the contents of file stdio.h, which is a 
‘header file for standard I/O‘. If mathematical functions such as cos and sin are 
used, we need the control line 

#include (math.h ) 


A program contains one or more functions. In the C language we have no 
sub-routines or procedures, just functions. Even the main program is a function, 
whose name is main. This is merely a matter of notation and terminology, since 
functions need not deliver a function value and can be invoked in exactly the same 
way as a procedure is called in other languages. This is shown in the following 
program, which reads the two-dimensional rectangular co-ordinates of two points P 
and О and computes the distance between these two points. 


/* This program computes the distance between 
two given points P and Q */ 


#include (math.h) 


main( ) 

( float xP, yP, xQ, уО; 
printf ("Give xP, yP, xQ, yQ:”); 
scanf ("96f %f of Pf”, &xP, &yP, &xQ, &yQ); 
print-distance (xP, yP, xQ, yQ); 


} 

print-distance (x1, y1, x2, y2) float x1, y1, x2, y2; 

{ float delta_x, delta-y, distance; 
delta x = x2 — x1; delta.y = y2 — yl; 
distance = sqrt (delta.x * delta x + delta.y * delta-y ); 
printf (^ Distance: Yof Nn", distance); 


The standard functions scanf, printf and sqrt will be discussed in Section A.9. 


A.5 ARRAYS AND POINTERS 


After the array declaration 
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the following five variables of type float are available: 
a[0], a[1], a[2], a[3], a[4] 


We can also write a[i], where subscript i can be any integer expression whose value 
will be neither negative nor greater than 4. Subscripts always count from 0. Between 
the brackets in the array declaration, only integer constants may appear. Instead of 
numbers, we often use named constants, as in the following example, which also 
shows that arrays can have more than one subscript: 


Y define NROWS 10 
#define NCOLUMNS 8 


int table/NROWS][NCOLUMNS]; 


for (1= 0; < NROWS; i++) 
for (j 0; j< NCOLUMNS; j++) 


„ЫД... 


If v is a variable, &v is the address of this variable, ог, in technical terms, &v is a 
pointer to v. If p is a pointer pointing to some object, that object is denoted by the 
expression *p. The following program shows how pointer variables can be declared 
and used: 


main ( ) 
{ inti, *p; 

i = 123; р = &i; xp = 789; 
} 


This is not exactly a practical program, but it shows that we can have access to a 
variable without using its name. Here p is a pointer to i, so *р is equivalent to i. 
This means that the variable i will finally have the value 789 instead of 123. 

The name of an array can be used as a pointer to its first element. Thus instead of 
&A[0] we can simply write A. If we have a pointer to an array element and we add 1 
to it we obtain a pointer to the next element (no matter how many bytes each 
element occupies). More precisely, instead of &A[i], we can simply write A + i. 

In their declarations arrays can be initialized if they have permanent memory 
space, that is, if they are either external or static. External variables are declared at 
the outermost level, that is, outside functions. The external arrays X and Y are 
initialized in the following program, which prints the value 5.75. 


float X[4] = (6.0, 6.0, 5.9, 6.1), 
Y[4] = (—0.25, 0.25, 0.0, 0.0); 


main( ) 
{ printf %f\n”, X[0] + Y[0]); 
} 
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Static variables have the keyword static in their declarations. In the following 
program the arrays X and Y are internal to the function main, but they are static, 
which means that they, too, have permanent memory space and can therefore be 
initialized. 

main( ) 
( static float X[4] — (6.0, 6.0, 5.9, 6.1), 
Y[4] = (—0.25, 0.25, 0.0, 0.0); 
printf ("?6f Nn", X[0]  Y[0]); 


A.6 FUNCTIONS 


In the C language a function may or may not have a function value. If it has, this 
value is delivered by means of a return-statement. As the name suggests, a 
return-statement causes an immediate return from the function. Thus the following 
two functions have the same effect. 


int f(x) float x; (if (x < 0) return —1; else return 1;} 
int f(x) float x; (if (x <0) return —1; return 1;} 


Incidentally, we can write a shorter version, using a conditional expression; 
int f(x) float x; (return x «0? —1:1;) 


Note how we express that a function has an integer function value and a float 
argument. The functions above can be called in the usual way, as, for example, in 


n=3*f(X*Y)+f(X); 


where X and Y are float variables. If an argument is declared float, we must not call 
the function with an integer argument. Thus if i is integer, f(i) is not a correct call of 
our function f. 

If a function does not deliver a function value it need not contain a return- 
statement. However, even then a return-statement can be used if we want an 
immediate return. We then simply omit the expression between return and the 
semicolon. 

In most cases it is only the numerical values of arguments that matter. However, 
sometimes we want a function that assigns values to variables, which are accessible 
through the arguments. This can be accomplished by letting these arguments be 
pointers to variables. The following function interchanges the values of two integer 
variables: 


interchange(p, q) int *p, *q; 
( intaux; 
} 


After the execution of 


і = 1; j = 2; interchange(&i, &j); 
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the values of i and j are 2 and 1, respectively. Note that the arguments &i and &j are 
pointers, which are denoted by p and q within the function. Two elements A[k] and 
A[m] of integer array A can be interchanged by this function in the following way: 


interchange(A +k, A + m); 


Remember that A + k is another notation for &A[k]. 

The above variable aux is said to have storage class automatic, which means that it 
has no permanent memory space. In contrast to this we have static and external 
variables, which are permanent. Automatic arrays cannot be initialized. A simple 
variable can always be initialized; if it is automatic, and declared in function F, it 
will have the initialization value each time F is entered. If it is static, however, it will 
have this value only the first time the function is entered. If the function is left and 
entered again later, a static variable will have its last value, whereas an automatic 
variable is undefined (if it is not initialized). 


A.7 STRUCTURES 


Several variables can be grouped together to a so-called structure. Suppose that we 
wish to store some information about points in two-dimensional space. For each 
point P this information consists of the rectangular co-ordinates x and y and a code 
to indicate whether the point lies inside a certain rectangle. Let us assume that there 
are the points P and Q. We can then declare 


struct (float x, y; int inside;) P, О; 


We can now use the components of structures P and Q in the same way as other 
variables, as, for example, in 


P-x=1.5;P-y=0.8; Р : inside = 1; 
Q:-x=2*P:x; 


Instead of the above declaration, we could have written either 


struct point {float x, y; int inside;} P, Q; 
or 
struct point {float x, y; int inside;}; 
struct point P, Q; 


These two versions have the advantage that the same structure type can be used 
later by simply writing struct point (as in the last declaration of P and Q), instead of 
struct {float x, y; int inside;}. Note that in this declaration we have to use the 
keyword struct. 

Should we wish to get rid of this obligation then we can use another facility, which 
in general enables us to give a type a new name. Here we can write 


typedef struct {float x, y; int inside;} POINT; 
Now POINT is a new name for our structure type, and we can proceed with 


POINT P, Q; 
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An array А of, say, 1000 points can be declared by 
struct (float x, y; int inside;) A[1000]; 
or, if the name point is defined as above, by 
struct point A[1000]; 
or, if typedef is used as above, by 
POINT A[1000]; 


Now each array-element A[i] is a structure consisting of the three components 
Ali] : x, Ali]: y, A[i] - inside. 

Structures are particularly useful in connection with dynamic memory allocation, 
since they can contain pointers to other structures. We can thus build lists, trees and 
so on. If S is a structure containing pointer field p then this pointer is $ - p and the 
object pointed to is *(S - p). For the latter expression a special notation is available, 
namely 


S—>p 


where the two characters — and > suggest an arrow, and form a single operator. 


A.8 DYNAMIC MEMORY ALLOCATION 
Suppose that we have declared 
char *p, *malloc( ), *realloc( ); 
Then p is a pointer to a character, and we can write 
p = malloc(n); 


where n is some positive integer expression, denoting a number of bytes. The effect 
is that, if possible, a consecutive piece of memory for n characters is allocated. If the 
required amount of memory is not available, p will have the value NULL. The latter 
is a special value for a pointer not pointing to a real object; it enables us to perform 
tests such as 


if (p == NULL) {... /* Insufficient memory */} 


The n bytes thus allocated are now available through pointer p. If, for example, the 
letter О is to be placed in the ith position (0x i € n — 1), we can write 


*(p +1) = 0"; 
The memory space thus allocated can be released by 
free(p); 


Suppose that a block of п bytes, pointed to by p, turns out to be too small, and 
that it is to be extended, the desired new size being N (>n). We can then write 


p = realloc(p, М); 
Now if p does not have the value NULL it will point to a block of N bytes; the first 
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n bytes of this block will have the same contents as those that p pointed to 
previously. 

We can also allocate memory space for other data types. Let us assume that we 
want a sequence of k integers. Since the standard function malloc must know how 
many bytes we need we are interested in the number of bytes occupied by an 
integer. This number is expressed in a machine-independent way by 


sizeof (int) 


Another complication concerns the pointer types that we now need. We need a 
pointer to an integer instead of a pointer to a character. However, malloc will yield 
a pointer to a character. Type conversion can be achieved by means of a 
cast-operator. Therefore we write 


int *p; 
p = (int *)malloc(k * sizeof (int)); 


Here the cast-operator (int *) says that the desired type is pointer to int. The jth 
integer (j = 0, 1,..., k — 1) of the sequence is denoted by *(p + j). 


. A.9 INPUT/OUTPUT 


In the C language there are no special constructs for input and output (I/O). 
Instead, we use a number of standard functions along with a predefined structure 
type, called FILE. Details about these are included in our program in a very 
convenient way, namely by writing 


#include (stdio.h ) 
The header file stdio.h contains a typedef declaration of the form 
typedef struct {...} FILE; 


Thus if we write 
FILE жўр; 


then variable fp has type pointer to FILE, and the compiler knows the details of the 
structure type FILE. Such a structure will actually be available through fp after a file 
is opened by 


fp = fopen(file-name, mode); 
where the arguments have the following meaning: 


file-name string containing the name of the file as it 
will appear in the directory; 

mode "r" for read, or 
"w" for write 


The counterpart of fopen is fclose; fopen connects a file to our program and fclose 
disconnects it. The following program shows how we can write something to the file 
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EXAMPLE. ЇЇ this file does not yet exist, it will be created. 


#include (stdio.h) 
main( ) 
( FILE *fp; 
int i; 
fp = fopen(’ EXAMPLE", "w"); 
for (i21;i«-4; it *) 
fprintf(fp, "i=%ld іжі = %2d\n”, i, i*i); 
fclose(fp); 
} 


After the execution of this program there will be the file EXAMPLE with the 
following contents: 


i=1 іжі = 1 
1-2 іжі= 4 
i23 ixi= 9 
i=4 ixi = 16 


Instead of fprintf we can use printf and omit the first argument, the output will 
then appear on the display of our workstation instead of on disk. The complete 
program would then read: 


main( ) 
{ int i; 
for (i=1;i<=4;i++) 
printf (^i = %1d ixi = %2d\n”, i, i*i); 


In fact printf(. . .) is equivalent to fprintf(stdout, ...), stdout being a file-pointer 
for standard output, declared in the header file stdio.h. The first argument of printf 
(the second argument of fprintf) is a format string containing pieces of text to 
appear literally in the output and format items, related to the remaining arguments. 
Here %14 and %2d are such format items; they specify that the values of i and i*i, 
are integer numbers, to be printed in one and two decimal positions, respectively. 
All other characters in the format string are printed literally. These characters 
include blank space and the newline-character, denoted by \n. In format items the 
letter d must be replaced with f if the associated data item is float instead of integer. 
Reading data is done similarly with 


fscanf(fp, format-string, . . .), 


if the data are in a file on disk, or 
scanf (format-string, . . .), 
if the data are to be typed in by the user. However, the remaining arguments must 


now be pointers, since fscanf and scanf must be able to assign values to variables 
through these arguments. For example, if a certain number must be given by the 
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user, we can write.: 
printf (" Give the number:"); scanf (° %d”, &number); 


(Omission of & in &number is a notorious error.) If x has type float and xx has type 
double, we can read their values as follows: 


scanf Pf 9olf", &x, &xx) 


The letter / in %/f is necessary because &xx has the type pointer to double. 

The function value returned by fscanf is the number of elements that have been 
read, so if this value is zero, nothing has been read. The latter property is often used 
to test whether the end of the file is reached. 

Analogous to stdout there is a standard file pointer stdin. The two calls scanf(. . .) 
and fscanf (stdin, . . .) are equivalent. 

A single character can be read and written in a more primitive way: 


ch — getc(fp) (input from disk) 
putc(ch, fp) (output to disk) 

ch — getchar( ) (input directly from user) 
putchar(ch) (output directly to user) 


Here getchar( ) is equivalent to getc(stdin) and putchar(ch) is equivalent to 
putc(ch, stdout). 

If we test the last character that has been read, and we decide that it is to be read 
again next time, we can put it back into the input stream by: 


ungetc(ch, fp) 


The functions printf, fprintf, scanf, fscanf perform so-called formatted I/O, as the 
final letter f in their names indicates. Formatted files have a line structure, and 
numbers are represented as sequences of characters. For files on disk that are only 
to be written and read by our own programs it is more efficient to use unformatted 
I/O. This means that the internal and the external representation of the data are 
identical, so numbers will probably be written as binary words of a fixed length. For 
unformatted I/O we use 


fread(bufptr, size, n, fp) 
fwrite(bufptr, size, n, fp) 


The arguments have the following types and meaning: 


bufptr pointer to char A pointer to a block of memory 
(sometimes called a buffer) 


size int Size in bytes of one data element 
n int Number of data elements in 
buffer 


fp pointer to FILE File-pointer 


These two functions deliver an integer as their function value, which is the number 
of data elements that has been read or written. As with fscanf, we can use this value 
to test whether we are trying to read beyond the end of the file, since then the read 
attempt will fail and the function value of fread will be zero. 
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А.10 MATHEMATICAL STANDARD FUNCTIONS 


Strictly speaking, the available set of predefined functions is not a part of the 
language. However, in practice it is convenient to have a list of such functions at 
hand. We have often used mathematical functions occurring in the following list, 
which shows function and argument types, along with a very brief indication about 
what the function computes. 


double cos(x) double x;  /*cos x */ 
double sin(x) double x; /* sin x */ 
double tan(x) double x; — /* tanx */ 
double log(x) double x; — /* lnx */ 
double sqrt(x) double x; — /* square root */ 
double floor(x) double x; /* floor (4.9) = 4.0 etc. */ 
double ceil(x) double x;  /* ceil (4.1) = 5.0 etc. */ 
int abs(i) int i; /* int absolute value */ 
double fabs(x) double x; /* double absolute value */ 
double acos(x) double x;  /* arccos x x/ 
double asin(x) double x; /* arcsin x ж/ 
double atan(x) double х; /* arctan x */ 
srand(seed) int seed ; /* initialize rand( ) */ 
int rand( ) /* random number generator */ 
long int time(p) /* time in seconds, since */ 


long int *p; /* 1 Jan. 1970 0.00h С.М.Т. */ 
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